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This  paper  presents  an  apertju  of  some  of  the  progress  in  the 
application  of  finite  element  methods  to  problems  of  continuum  and 
structural  mechanics.  Since  space  was  limited,  we  could  not  repro¬ 
duce  the  theory  and  applications  of  an  extensive  class  of  new  equil¬ 
ibrium  elements,  and  had  to  restrict  the  review  to  pure  kinematic 
models  {matrix  displacement  method).  Moreover,  in  this  latter  field 
we  can  only  give  a  limited  survey  of  the  work  carried  out  at  the 
Institute  of  Statics  and  Dynamics  and  the  Department  of  Aeronautical 
Structures  at  Imperial  College.  With  few  exceptions,  we  have  ignored 
the  basic  theory  and  have  reproduced  only  some  characteristic  results. 
Iterative  methods  of  solutions  of  problems  with  a  large  number  of 
unknowns  using  the  conjugate  gradient  method  are  discussed  in 
References  20  and  29. 

Broadly,  the  research  presented  here  covers  the  following  fields: 
a)  Application  of  more  refined  elements  to  three-dimensional  stress 
and  strain  fields;  reference  to  the  corresponding  two-dimensional 
developments  are  not  made,  b)  Application  of  two  highly  efficient 
classes  of  fully  compatible  plate  elements,  c)  A  new  class  of  elements 
for  shells  of  arbitrary  geometry,  d)  A  generalization  and  refinement 
of  the  analysis  of  elastro-plastic  deformations  of  Reference  1,  in¬ 
cluding  a  new  iterative  method;  (References  23  and  24).  Here  too,  we 
had  to  omit  a  discussion  of  the  kinematic  hardening  effect  (Refer¬ 
ence  25).  e)  Some  techniques  for  the  analysis  of  large  displacements 
and  strains. 

We  cannot  enter  the  far  reaching  evolution  of  the  structural 
language  ASKA  (Reference  27),  which  bears  but  little  connection  with 
the  ASKA  packet  discussed  in  Reference  1. 


♦Professor,  Director 
**Group  Leader 
***Deputy  Group  Leader 
♦♦♦♦Scientific  Collaborators 
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SECTION  I 

THREE-DIMENSIONAL  ANALYSIS 

1.  GENERAL  CONSIDERATIONS 

The  analysis  of  two  and  three-dimensional  stress  and  strain  fields  was  discussed  at 
considerable  length  in  a  paper  at  the  First  Dayton  Conference  (Reference  1).  In  the  intervening 
three  years  significant  advances  have  been  achieved  in  evolving  new  elements  of  greater 
generality  and  refinement. 

In  this  connection  two  associated  aspects  received  particular  attention.  The  first  concerns 
a  more  faithful  reproduction  of  the  actual  geometry  of  a  structure  than  is  feasible  with  the 
original  set  of  elements,  which  demands,  in  general,  a  polygonization  or  polyhedrization  of 
the  system.  For  some  time  it  has  been  evident  that  by  ignoring  the  actual  curvature  involves, 
in  certain  cases,  a  marked  loss  of  accuracy.  The  second  extension,  closely  allied  to  the  first 
one  realized  in  the  theory,  is  the  introduction  of  polynomials  of  higher  order  for  the  dis¬ 
placement  fields.  Moreover,  while  in  the  past,  pure  kinematic  compatibility  was  aimed  at, 
some  of  the  new  elements  allow  also  for  continuity  in  the  first  derivatives  at  certain  nodal 
points.  This  means  that  stresses  and  rotations  may  be  taken  as  continuous  at  these  nodal 
points  if  this  is  physically  admissible.  No  doubt  this  yields  a  more  accurate  stress  field, 
which  is  also  easier  to  interpret.  However,  the  incorporation  of  more  degrees  of  freedom 
than  are  strictly  necessary  from  the  variational  formulation  of  the  problem,  involves  a 
so-called  overcompatibility  with  a  consequent  slight  increase  in  the  stiffness  of  the  element. 
Some  of  the  relevant  progress  in  this  field  is  discussed  briefly  in  References  3  to  6. 

Broadly  speaking,  it  may  be  stated  that  the  new  elements  are  based  on  the  application  of 
a  number  of  distinct  interpolation  techniques  for  the  direct  generation  of  the  so-called  ‘unit’ 
displacement  functions.  The  latter  leads  to  the  immediate  formation  of  the  elastic  stiffness 

and  may  be  also  used  in  a  simple  manner  to  obtain  the  geometrical  stiffness  .  For 
example,  we  have  the  wide  class  of  LUMINA  elements,  which  evolve  from  Lagrange’s  inter¬ 
polation  technique.  These  elements  have  a  three-dimensional  hexahedronal  form,  whose  six 
faces  are  curved  in  contrast  to  what  is  commonly  called  a  hexahedron  element  (Figure  1). 
The  geometry  of  the  element  is  fixed  through  a  unique  one  to  one  mapping  of  the  physical 
space  0  x  x  x  onto  a  regular  cube  I  £  I  <  I  in  the  0  £  £  £  space  by  Lagrangian  interweaving 
interpolation  functions  through  a  selected  set  of  nodal  points  (Figure  2).  A  systematic  and 
Consistent  arrangement  of  nodal  points  Is  achieved  by  prescribing  (1  +  1),  (m  +  1),  (n+1)  points 
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(  n  +1 )  Nodol  Points 


Figure  1.  LUMINA  -  A  General  Hexahedron  Element 


last 


/  *  ( I  +  1)(  ititt]  (n+l) 

Figure  2.  LUMINA  -  Basic  Cube  Nodal  Point  Numbering 
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in  the  three  curvilinear  directions  of  the  actual  element.  Thus  we  note  that  there  are  in  all 
r  =  (1+1)  (m+1)  nodal  points.  As  far  as  the  selection  of  the  fixed  set  of  nodal  points  is  con¬ 
cerned,  the  reader  should  consult  Reference  3.  The  order  of  the  Lagrangian  interpolation  is 
clearly  1,  m,  n  respectively  in  the  basic  mapping  space.  It  is  evident  that  these  functions  are 
not  complete  in  the  £  space.  When  1  =  m  =  n,  the  elements  in  question  are  denoted  within 
ASKA  by  HEXE  r.  For  an  interpolation  function  of  the  second  order  we  note  that  the  hexa- 
hedronal  element  has  27  nodal  points  and  81  unknowns.  For  a  third  order  interpolation 
function  in  all  directions  we  have  HEXE  64  with  192  unknowns.  These  elements  are  used 
extensively  in  the  analysis  of  turbine  blades,  arch  dams  and  other  three-dimensional  con¬ 
figurations  where  the  curvature  is  a  decisive  factor.  It  is  of  interest  that  HEXE  r  has  a 
complete  displacement  field  of  up  to  the  r  order  in  the  £  -space. 

Another  three-dimensional  interpolation  element  in  ASKA  is  called  HERMES  8  is  dis¬ 
cussed  in  Reference  4.  Just  as  in  the  LUMINA  set,  the  HERMES  element  represents  a  general 
hexahedronal  element  with  curved  faces,  and  has  proved  a  most  useful  component  block  for 
three-dimensional  analysis  of  complex  bodies.  The  cardinal  idea  underlying  the  HERMES  de¬ 
velopment  aims  at  combining  the  advantages  of  the  Lagrangian  and  Hermitian  interpolation 
techniques.  These  maybe  summarized  as  follows: 

Lagrangian  Interpolation  Technique.  -  The  geometry  of  the  element  is  described  most  easily 
and  systematically  by  an  essential  set  of  nodal  points,  no  derivatives  or  normal  vectors  to  the 
surfaces  are  required. 

Hermitian  Interpolation  Technique.  -  The  derivatives  -  and  hence  strain  and  rotations  -  enter 
directly  as  independent  parameters  defining  the  stiffness  and  other  relevant  properties  of  the 
element.  This,  if  physically  admissible,  allows  one,  as  stated  above,  to  stipulate  continuity  of 
strains  and  stresses. 

It  should  be  expected  that,  for  the  same  number  of  freedoms,  the  resulting  strains  and 
stresses  generally  will  be  more  accurate  when  using  HERMES  elements  than  in  the  straight¬ 
forward  application  of  LUMINA  elements,  where  the  strains  (and  hence  stresses)  are  deduced 
by  way  of  the  derivatives  of  modal  Lagrangian  functions. 

HERMES  8,  which  is  the  simplest  three-dimensional  Hermitian  interpolation  element, 
has  all  its  freedoms  assigned  to  the  eight  corner  points,  hence  any  intermediate  nodal  points 
necessary  for  the  geometrical  shaping  of  the  element,  as  in  LUMINA,  do  not  have  any 
freedoms  and  are  dummy  points  as  far  as  the  kinematic  analysis  is  concerned. 
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B  Is  evident  that  HERMES  8  uses  a  partial  set  of  third  order  Hermltian  polynomials.  For 
all  the  usefulness  and  considerable  success  of  the  LUMINA  and  HERMES  elements,  there  is 
little  doubt  that,  for  arbitrary  stress  fields,  it  is  desirable  to  apply  elements  with  complete 
polynomials  for  the  displacement  fields.  Following  Reference  30  this  is  in  practical  terms 
only  possible  with  triangular  and  tetrahedron  elements  for  two-and  three-dimensional  stress 
fields  respectively.  In  the  past  such  elements  (see  e.g.  TRIM  6  and  TET  10  of  Reference  1) 
have  been  evolved  allowing  up  to  a  linear  variation  of  strains,  i.e.  a  quadratic  variation  of 
displacements.  In  general  these  proved  satisfactory,  but  the  necessity  arose  to  cater  for  a 
greater  sophistication  in  the  specification  of  the  displacements.  This  is  achieved  by  the 
TET  20  and  the  TEA  8  elements  presented  in  Reference  5.  B  has  been  shown  in  Reference  30 
and  stated  in  Reference  2  that  third  order  complete  polynomials  yielding  a  second  order  or 
parabolic  strain  distribution  fit  Into  a  tetrahedron  element  with  20  nodal  points  and  60  degrees 
of  freedom,  denoted  in  ASKA  as  TET  20  Figure  3.  The  nodal  points  of  this  element  include,  in 
addition  to  the  four  vertices,  two  points  on  each  edge  trisecting  the  length  -  a  total  of  twelve 
points  -  and  four  points  at  the  centroids  of  the  faces.  The  geometrical  choice  of  the  edge  nodal 
points  and  centroids  is  dictated  by  symmetry  and  algebraic  convenience.  This  new  element 
satisfies,  of  course,  kinematic  compatibility  of  the  displacements  everywhere. 


An  interesting  variant  of  the  use  of  third  order  polynomials  presents  iself,  if  we  stipulate 
that  kinematic  compatibility  at  the  four  vertices  should  be  extended  to  all  first  derivatives 
of  the  displacements.  The  associated  degrees  of  freedom  are  clearly  4  x  12  =  48.  The  nodal 
points  of  the  edges  are  then  redundant,  but  to  satisfy  kinematic  compatibility  of  the  dis¬ 
placements  over  the  faces  of  the  tetrahedron  it  is  necessary  to  retain  the  four  ‘standard’ 
nodal  points  at  the  centroids  of  the  faces  of  the  TET  20,  as  shown  in  Figure  4.  The  total 
degrees  of  freedom  of  such  an  element  -  called  TEA  8  in  ASKA  -  Is  clearly  60,  as  for  the 
TET  20. 


It  is  Intriguing  to  observe  that  application  of  the  modal  technique  of  Reference  1  leads 
to  a  so-called  basic  or  modal  stiffness  k  of  the  order  (60  x  60),  from  which  the  stiffness 
of  TET  20  and  TEA  8  elements  may  be  derived  from  two  alternative  congruent  transformations. 

Essentially,  the  technique  of  Reference  5  Is  based  on  the  modal  idea  of  Reference  1. 
However,  a  more  general  approach  may  be  achieved  by  a  direct  construction  of  the  dis¬ 
placement  fields  on  the  basis  of  a  distinct  and  complete  interpolation  technique. 
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Consider  a  general  TET  or  TRIM  element  corresponding  to  a  complete  polynomial  of 
th 

nr  order.  Such  a  construction  is  always  feasible,  as  long  as  the  necessary  number  of  internal 
nodal  points  f.  are  allowed  for.  We  note  that 


fj  =  ~-(m+l)(m+2) 

for  TRIM  elements 

(  II 

f:  =  ~  (m+l)(m  +  2)(m+  3) 

for  TET  elements 

(2) 

The  modal  functions  o> ,  corresponding  to  unit  displacement  at  one  nodal  point  and  zero  at  all 
others,  may  now  be  determined  directly  without  first  establishing  a  transformation  matrix.  It 
is  only  necessary  to  assume  a  regular  array  of  nodal  points  and  apply  the  complete  inter¬ 
polation  scheme  in  £  shown  in  Figure  5.  This  procedure  has  been  applied  within  ASKA  to 
TRIM  15,  which  corresponds  to  m  =  4.  Similar  techniques  may  be  developed  for  the  corre¬ 
sponding  neo-dimens ional  element  groups  TRIA  and  TEA  of  the  order  m  >  3.  For  example, 
there  is  the  element  TRIA  6  for  m  =  4,  which  has  three  ‘special’  nodal  points  at  the  vertices 
and  three  ‘simpler’  ones  at  the  mid-points,  where  the  displacements  and  their  derivatives 
normal  to  edges  are  prescribed;  (15  modal  functions),  cubic  variation  of  strain. 

The  procedure  may  also  be  amplified  for  m>  2  to  TET  and  TRIM  elements  with  curved 
boundaries  (edges  and  faces).  It  is  only  necessary  to  use  the  same  interpolation  method  of 
Figure  5,  as  a  one  to  one  mapping  device,  to  construct  the  actual  curvilinear  triangular  and 
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Figure  5.  General  TET  and  TRIM  Elements  -  Complete  Interpolation  Scheme 
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tetrahedronal  elements  from  the  corresponding  basic  configurations  with  straight  edges  and 
faces.  Such  elements  are  denoted  within  ASKA  by  TRIC  N  and  TEC  N  respectively,  N  rep¬ 
resenting  the  number  of  nodal  points.  Moreover,  the  definition  of  a  geometrical  stiffness 
proves  simple  now  for  these  elements  and  the  technique  may  be  extended  to  the  large  dis¬ 
placement  analysis  of  two-  and  three-dimensional  systems  (Reference  6). 

2.  APPLICATIONS 

a.  Analysis  of  an  Archdam  Using  HEXE  64  Elements 

The  geometry  of  the  archdam  analyzed  here  is  shown  in  Figure  6.  The  curved  surfaces 

of  this  system  can  be  approximated  best  by  interpolation  elements,  such  as  LUMINA  or 

HERMES.  In  the  present  case  the  idealization  uses,  one  layer  of  18  curved  HEXE  64  elements, 

the  grid  of  which,  projected  onto  the  reference  cylinder,  is  given  in  Figure  6.  The  material 

6  2 

constants  E  and  v  are  taken  as  0.20  x  10  kp/cm  and  0, 15  respectively.  A  kinematically 
consistent  idealization  of  the  water  pressure  loading  is  used  in  accordance  with  the  deformation 
modes.  The  base  of  the  dam  is  assumed  to  be  built  in  rigidly  with  zero  displacements  at  the 
nodal  points.  This  extreme  discretization  was  adopted  at  the  specifc  request  of  the  Civil 
Engineering  Research  Council  at  the  Institution  of  Civil  Engineers,  Great  Britain,  for  whom 
this  analysis  was  carried  out. 

Some  typical  results  of  the  computations  for  the  air  and  water  face  stresses  at  the  plane 
of  symmetry  are  shown  in  Figures  7  and  8.  On  examining,  for  instance,  the  radial  stresses 
at  the  water  face,  these  will  be  seen  to  correspond  very  closely  indeed  to  the  applied  water 
pressure,  except  near  the  base  of  the  dam,  where  the  effect  of  the  rigid  building-in  is  clearly 
apparent.  These  stress  peaks  can  be  smoothed  out  using  a  foundation  model  with  six  additional 
HEXE’s  representing  the  rock  bed.  For  details  of  this  calculation  see  Reference  7.  A  corre¬ 
sponding  TET  10  analysis  of  the  archdam  described  here  is  found  in  Reference  8,  and  shows 
an  excellent  agreement  with  the  current  results. 

b.  Analysis  of  a  Heated  Turbine  Vane  Using  HERMES  8  Elements 

The  HERMES  8  element  has  been  applied  to  the  analysis  of  the  cylindrical  vane  with  two 
cooling  channels  enclosed  by  solid  end  parts,  as  shown  in  Figure  9.  This  structure  is  of  con¬ 
siderable  complexity  not  only  as  far  as  the  geometry  is  concerned  but  also  with  respect  to 
the  loading  conditions.  The  temperature  varies  both  in  cross-sectional  and  axial  directions, 
but  is  symmetrical  about  the  middle  cross-section,  as  is  the  geometry.  Therefore  we  may 
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Figure  7.  Air  Face  Stresses  at  Plane  of  Symmetry  Figure  8.  Water  Face  Stresses  at  Plane  of  Symmetry 
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only  consider  half  of  the  system  suppressing  at  the  ‘root’  all  displacements  orthogonal  to 
the  plane  of  symmetry.  Due  to  the  change  of  temperature  within  a  large  range,  it  is  also 
necessary  to  allow  for  variation  of  the  material  properties  as  Young’s  modulus,  Poisson’s 
ratio  and  thermal  coefficients  of  expansion. 

The  vane  is  subdivided  into  a  suitable  number  of  hexahedra  with  plane  surfaces  and  side 
length  ratios  of  approximately  1:3:6.  The  typical  middle  cross  section  of  the  three- 
dimensional  grid  is  shown  in  Figures  10  and  11.  The  question  is  why  not  a  more  general 
HERMES  element,  with  the  desirable  feature  of  curved  surfaces?  The  reason  is  that  the  same 
vane  had  been  calculated  previously  with  six  TET  10  elements  (Reference  1)  in  place  of  one 
hexahedron  (Reference  9),  In  order  to  check  the  agreement  between  both  computations  we 
aimed  in  changing  the  geometry  as  little  as  possible.  To  simplify  the  HERMES  8  idealization 
a  slight  modification  had,  however,  to  be  introduced  near  the  trailing  edge  (Figure  11).  In  the 
TET  10  analysis  the  material  properties  are  taken  as  constant  within  each  tetrahedron  using 
a  mean  temperature,  but  the  thermal  strains  were  allowed  to  vary  linearly.  Within  a  HERMES  8 
element  Young’s  modulus  and  Poisson’s  ratio  as  well  as  the  thermal  strains  are  inter¬ 
polated  using  the  Lagrangian  polynomial  set  of  HEXE  27. 

It  is  most  interesting  to  compare  the  number  of  elements  and  unknowns  involved  in  the 
TET  10  and  HERMES  8  discretizations 


Number  of 

TET  10 

HERMES  8 

Elements 

906 

155 

Nodal  points 

1770 

304  (1825) 

Unknowns 

5192 

3531 

Due  to  the  adopted  HEXE  27  interpolation  scheme  the  HERMES  8  idealization  requires 
the  temperatures,  material  data,  and  initial  strains  at  1825  nodal  points.  However,  kinematic 
freedoms  are  only  assigned  to  304  nodal  points.  As  the  96  degrees  of  freedom  of  a  HERMES  8 
element  are  concentrated  in  8  nodal  points,  which  are  common  for  many  adjacent  elements, 
the  number  of  unknowns  is  remarkably  reduced.  For  lack  of  space  we  restrict  our  presentation 
to  some  characteristic  results  Figures  12  and  13  show  a  comparison  of  the  spanwise  dis¬ 
placement  w  as  calculated  with  TET  10  and  HERMES  8  elements.  The  distribution  of  the 
greatest  component  stress  crzz  is  plotted  in  Figures  14  and  15  for  two  typical  cross  sections. 
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Figure  10.  Basic  Grid  of  Elements 


Figure  11.  Modification  of  Grid  for  HERMES  8  Elements 
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Figure  13.  Comparison  of  Displacement  w  -  TET  10  and  HERMES  8  Analyses 
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Figure  15.  Typical  Distribution  of  Component  Stress  <rzz  i 
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SECTION  n 

ANALYSIS  OF  THIN  PLATES 

1.  GENERAL  CONSIDERATIONS 

The  analytical  construction  of  displacement  functions  for  plate  elements  subject  to  bending 
was  discussed  at  some  length  in  Reference  2.  Attention  was  drawn  to  both  triangular  and 
quadrilateral  elements  available  in  the  ASKA  system,  which  satisfy  either  all  kinematic 
compatibility  conditions  or  are  deficient  in  the  continuity  of  the  gradient  normal  to  the  edge. 
The  most  rudimentary  elements  of  the  triangular  set  consist  of  the  so-called  fully  compatible 
TRIB  3C,  and  the  TRIB  3  (References  1,  and  18),  which  satisfies  the  slope  condition  only  at 
the  vertices.  The  TRIB  3C  corresponds  in  principle,  but  not  in  detail,  to  the  element  evolved 
by  Bazeley  et  al  (Reference  15),  but  allows,  however,  for  linear  taper.  The  comments  on  the 
pros  and  cons  of  these  elements  were  summarized  in  Reference  2.  Note  also  the  interesting 
construction  given  by  Clough  and  Tocher  (Reference  16)  and  the  elegant  solution  to  the  matching 
of  two  triangular  elements  developed  by  de  Veubeke  (Reference  17). 

It  was  also  pointed  out  in  Reference  2  that  a  more  systematic  examination  of  the  funda¬ 
mental  problems  arising  in  triangular  elements  under  bending  leads  to  the  conception  of  a 
number  of  distinct  families  for  such  elements,  which  are  fully  compatible  and  consistent  in 
their  convergence.  While,  in  principle,  more  sophisticated  than  TRIB  3C,  they  are  essentially 
simpler  and  more  accurate  in  their  final  formulation.  We  discuss  briefly  one  particular 
family  known  in  the  ASKA  system  as  the  TUBA  set. 

The  TUBA  family  is  based  on  complete  polynomial  functions  for  the  deflection  w  of  the 

order  mi  5.  The  selection  of  nodal  freedoms,  which  include  at  the  vertices  not  only  w,  and 

the  first  derivatives  w,  ,  w,  but  also  all  the  second  derivatives  w,  ,  w,  ,  w,  ,  is 

x  y  xx  xy  yy 

characteristic.  Additional  ‘simpler’  nodal  points  are  strategically  placed  on  the  boundary 
to  ensure  complete  compatibility  for  w  (order  m)  and  the  slope  normal  to  the  edges  w,n 
(order(m-l)).  If  necessary,  further  freedoms  f^  are  assigned  to  interior  points  to  define 
uniquely  the  s  constants 

$  =  -g-  {m  +  I)  (  m  +2  )  (3  ) 

of  the  polynomial  or  order  m.  It  is  easy  to  confirm  that  such  a  scheme 

fj  =-£  (m  -4  )  (  m-5  )  (4) 
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Thus,  for  a  polynomial  of  order  m=  5,  there  are  no  interior  points,  the  corresponding  element 
being  called  TUBA  6,  since,  in  addition  to  the  three  vertices,  three  nodal  points  have  to  be 
placed  in  the  boundary  (Figure  16).  For  m  =  6  we  require,  in  addition  to  12  nodal  points  on 
the  edge,  one  interior  point,  the  element  being  denoted  by  TUBA  13  (Figure  17).  The  highest 
order  of  TUBA  element  available  in  ASKA  uses  m  =  7,  for  which  f.  =  3,  the  connotations  of 
the  element  being  TUBA  15  (Figure  18).  No  practical  necessity  has  arisen  until  now  to  expand 
the  TUBA  group  with  higher  order  functions,  although  the  procedure  is  straightforward.  The 
number  of  kinematic  freedoms  for  TUBA  5,  13,  15  is  seen  to  be  21,  28,  36  respectively. 

One  may,  of  course,  criticize  the  basic  philosophy  of  these  elements  for  introducing,  as 

freedoms,  the  three  second  order  derivatives  at  the  vertices,  which  do  not  appear,  in  general, 

in  a  variational  formulation  of  the  problem,  and  lead  to  the  aforementioned  overconformity. 

On  the  other  hand,  the  prescription  of  the  derivatives  w,  ,  w,  ,  w,  at  the  vertices  ensures 

ay  jj 

the  continuity  of  moments  or  stresses  there,  if  this  is  admissible,  which  proves  valuable  from 
the  practical  point  of  view,  since  it  yields  accurate  estimates  of  these  measures.  Moreover, 
the  selected  idealization  fits  neatly  into  the  general  concept  of  the  ASKA  system.  (Ref¬ 
erences  12  and  13)  describe  the  main  outline  of  the  theory  for  the  TUBA  6  and  TUBA  15 
elements.  The  TUBA  set  has  proved  highly  efficient  in  the  analysis  of  plates  of  arbitrary 
shapes  under  static  and  dynamic  loading.  This  is  illustrated  by  the  examples  given. 

However,  for  the  specific  geometry  of  parallelogram  plates,  it  may  be  advantageous  to 
use  different  complete  sets  in  oblique  coordinates.  To  this  purpose  we  have  developed  PUBA  4, 
using  third  order  Hermitian  polynomials  and  allowing  a  complete  displacement  function  of  w 
of  the  third  order.  As  freedoms  we  select  the  w,  w,x>  w,y,  w>xy  at  the  four  nodal  points 
(16  freedoms).  Extension  to  higher  order  Hermitian  interpolation  polynomials  is  evident.  The 
complete  theory  of  PUBA  elements  is  developed  in  Reference  14;  they  may  be  considered  as 
a  very  broad  generalization  of  the  rectangular  elements  evolved  by  Bogner,  Fox  and  Schmit 
(Reference  21).  They  allow,  in  addition  to  arbitrary  obliquity,  also  for  arbitrary  anisotropy 
of  the  material.  Of  importance  is  that  for  both  sets  of  elements  the  geometrical  stiffness 
for  given  in-plane  membrane  loads  has  been  set  up.  This  allows  the  investigation  of  instability 
phenomena  in  plates. 

2.  APPLICATIONS 

a.  The  TUBA  6  Element 

In  order  to  check  the  .ccuracy  of  the  TUBA  6  element  we  use  the  classical  test  example 
of  a  simply  supported  square  plate  homogeneous  and  isotropic,  v  -  0.3).  Due  to  the  symmetry 
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a)  Natural  Freedoms  pTM  ..  _  .  . 

rTN  b)  Cartesian  Freedoms  pj 

Figure  16.  TUBA  6  Element 


•Natural  First  Derivative  dw/dt,- 


— >  Natural  Derivative  d  <u/dv  =  w 
resp  slope  =  dw/dn  =w 

''  «  A 

— >Nafurol  Second  Derivative  dtw/dCi 

j 

a)  Natural  Freedoms  pTN 


w , .  . 

i ) 


Prescri  be  d 
•  Deflection  w 
o  First  Derivative  (  s  ) 

(*)D  ef  I  e  ct  ion,  First,  Second  Derivatives 


b)  Cartesian  Freedoms 


Pj 


Figure  17.  TUBA  13  Element 
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of  the  system  we  may  only  consider  one  quarter  of  the  plate;  strictly,  this  could  be  reduced 
to  one  eighth.  It  is  first  necessary  to  introduce  the  kinematic  boundary  conditions,  e.g.  to 
suppress  at  a  support  parallel  to  a  y-axis  the  freedoms  w,  w,  ,  w,  .  Then,  in  order  to 

y  yy 

satisfy  the  dynamic  boundary  conditions,  w>xx  is  nullified  at  such  a  location.  Alternatively 
and  equivalently,  we  may  introduce  zero  moment  m  as  an  external  loading.  Two  loading 
cases  are  analyzed,  a  concentrated  load  at  the  center  and  a  uniform  pressure  represented  by 
kinematically  equivalent  loads.  The  excellent  agreement  of  the  finite  element  results  with 
the  exact  solution  is  apparent  in  Figure  4.  Some  of  the  trigonometric  series  of  the  analytical 
theory  have  been  recalculated  to  a  higher  accuracy  than  found  in  the  standard  tables  in  order 
to  have  a  more  reliable  yardstick  for  our  computations.  The  only  remarkable  deviation  is 
found  in  the  case  of  the  concentrated  load  for  the  bending  moment  at  a  point  B  of  Figure  19, 
which  lies  rather  close  to  the  singularity.  There  is  little  doubt  that  TUBA  6  represents  the  best 
plate  element  available  to  date.  We  regret  that  lack  of  space  does  not  allow  any  comparison 
with  TUBA  13  and  TUBA  15  results. 

b.  The  PUBA  4  Element 

A  convergence  test  similar  to  that  of  the  last  section  has  been  undertaken  applying  PUBA  4 
elements  for  the  idealization  of  a  skew  plate.  For  the  deflection  analysis  we  consider  a  plate 
with  a  skew  angle  of  60*;  note  that  the  skew  angle  for  a  rectangular  plate  is  here  90*.  The 
discretization  of  the  plate  proceeds  with  4,  8,  12  and  16  elements  per  side.  Kinematically 
equivalent  loads  representing  a  uniform  pressure  are  applied.  Following  Figure  20  the  per¬ 
centage  error  in  the  maximum  deflection  with  respect  to  the  analytical  values  given  by 
Morley  (Reference  37)  decreases  from  +  1.65%  to  0.4%,  which  is  quite  remarkable.  We  also 
analyze  the  buckling  of  skew  plates  subject  to  an  oblique  longitudinal  stress.  In  this  case  it  is 
senseless  to  draw  curves  confirming  the  convergence,  because  with  four  elements  per  side 
a  numerical  accuracy  of  three  decimal  digits  is  achieved.  We  therefore  limit  the  presentation 
in  Figure  21  to  the  variation  of  the  critical  stress  versus  the  skew  angle. 
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Figure  18.  TUBA  15  Element 
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Figure  20.  PUB  A  4  Element,  Convergence  of  Maximal  Deflection  of  Simply  Supported 
Skew  Plate 


/’ 


Figure  21.  PUBA  4  Element,  Critical  Stress  for  a  Built-in  Plate  as  Depending  on  the 
Skew  Angle 
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SECTION  III 

THE  SHEBA  SHELL  ELEMENT  FOR  THE 
MATRIX  DISPLACEMENT  METHOD 

1.  GENERAL  CONSIDERATIONS 

Considerable  effort  has  been  extended  over  the  last  years  in  adapting  the  matrix  dis¬ 
placement  method  to  the  specific  problems  of  thin  shells  under  membrane  and  bending  action 
and  developing  suitable  elements  varying  in  sophistication.  Some  of  the  difficulties  arising  in 
the  process  of  idealization  were  reviewed  in  Reference  2.  For  example,  simple  considerations 
show  that  a  representation  of  a  shell  by  polyhedron  surfaces  may  lead  to  serious  errors, 
especially  in  the  presence  of  pronounced  bending  and  so-called  boundary  layer  effects.  For 
this  and  other  reasons  it  appears  imperative  to  allow  for  the  curvature  of  the  shell.  It  is, 
however,  evident,  even  from  a  cursory  perusal  of  the  literature,  that  no  reliable  procedure 
exists  for  shells  of  arbitrary  geometry  and  especially  systems  with  negative  Gaussian  curva¬ 
ture  which  are  of  increasing  importance  for  cooling  aggregates  and  similar  systems.  More¬ 
over,  much  of  the  work  developed  to  date  is  deficient  as  far  as  some  basic  requirements  of 
the  matrix  displacement  method  are  concerned.  For  example,  It  is  well  known  (Reference  1) 
that  one  fundamental  prerequisite  for  establishing  reliable  stiffness  matrices  is  the  existence 
of  true  rigid-body  motions,  giving  rise  to  zero  strains.  Without  this  condition  the  necessary 
equilibrium  of  an  element  may  not  be  satisfied,  which  may  lead  to  serious  errors.  It  is  now 
well  established  with  the  classical  work  of  Reissner  (Reference  35),  Koiter  (Reference  33), 
and  Sanders  (References  32  and  34),  that  Love’s  formulation  of  the  thin  shell  theory  is  deficient 
since  it  gives  rise  to  strains  under  rigid-body  rotation  unless  the  shell  is  spherical.  The 
effect  is  also  nonexistent  in  shells  of  revolution  under  symmetrical  loading.  The  error  is,  of 
course,  not  serious  as  stressed  by  Reissner,  Koiter  and  Sanders,  as  long  as  the  shell  is  not 
subjected  to  pronounced  twisting  action.  Nevertheless,  this  conclusion  applies  strictly  only  to 
an  analysis  of  a  shell  as  a  single  system  but  loses  much  of  Its  validity  when  we  have  to  con¬ 
sider  a  finite  element  idealization.  Here  equilibrium  has  to  be  satisfied  exactly  for  each 
element.  Moreover,  any  extension  into  the  realm  of  large  displacements  becomes  of  question¬ 
able  reliability  if  this  criterior  is  not  satisfied.  Any  simplifications  of  the  basic  theory  as 
expressed,  for  example,  by  the  shallow  shell  formulation  of  Marguerre- Vlas  sow  have  also 
to  be  viewed  critically,  when  applied  to  finite  elements.  This  is  particularly  true  when  con¬ 
sidering  kinematic  compatibility  between  adjoining  elements. 

With  this  background  in  mind  and  with  the  necessity  to  provide  an  element  applicable  to 
any  geometry  and  shell  configuration  we  observe  once  more  that  the  most  useful  element 
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must  be  a  triangular  one,  with  curvilinear  boundaries  and  varying  curvature.  In  fact,  this 
should  represent  a  very  broad  generalization  of  the  TUBA  set  of  Reference  12.  It  is  shown 
in  the  present  article  that  interelement  compatibility  can  be  achieved,  if  the  three  global 
cartesian  displacements  are  connected  to  the  kinematic  freedoms  at  the  nodal  points  by  the 
same  transformation  or  interpolation  procedure  as  used  for  the  deflection  w  in  the  TUBA 
family  (Reference  12).  Moreover,  in  order  to  ensure  zero  strains  under  rigid-body  motion 
it  is  necessary  to  develop  a  mapping  of  the  surface  in  global  cartesian  coordinates  based 
on  the  corresponding  set  of  nodal  geometrical  data  and  obeying  the  same  rule  of  transfor¬ 
mation  as  the  displacements.  The  transformation  matrices  are  expressed  most  conveniently 
in  terms  of  the  homogeneous  triangular  coordinates  referred  to  the  so-called  basic  plane 
triangle  formed  by  the  vertices  of  the  element. 

One  aspect  that  arises  immediately  in  the  initiation  of  such  an  investigation  is  the  question 
whether  it  is  not  most  convenient  to  operate  with  displacements  in  the  directions  of  the 
principal  curvatures  and  the  normal  to  the  surface,  as  in  the  standard  shell  theory.  However, 
since  the  mapping  of  the  geometry  and  the  identical  one  for  the  displacements  is  expressed 
for  the  cartesian  components  in  terms  of  the  natural  homogeneous  triangular  coordinates  of 
the  basic  triangle.  It  is  found  that  the  transformations  to  the  directions  of  principal  curvature 
become  extremely  complicated  and  inconvenient,  especially  as  far  as  the  second  order 
derivatives  are  concerned.  It  soon  becomes  evident  that  it  is  best  to  keep  within  the  ‘natural’ 
concept  originated  in  Reference  1  and  define  both  membrane  strains  and  curvatures  in  the 
directions  of  the  so-called  a  ,  /3,  y  curves  corresponding  to  C, ,  =  const,  respectively. 

The  theory  may  now  be  evolved  in  a  most  elegant  manner  and  leads,  without  any  artifices,  to 
a  consistent  formulation  yielding  true  zero  strains  under  any  rigid-body  motion.  This  is  truly 
remarkable  in  view  of  the  basic  simplicity  of  the  analysis.  Naturally,  our  considerations 
simplify  considerably  for  shell  elements  of  specific  geometry  such  as  circular  cylindrical 
shells,  shells  of  revolution  and  spherical  shells. 

The  general  theory  is,  however,  necessary  from  another  point  of  view.  Even  for  a  shell 
of  simplified  geometry  generality  is  required  once  we  have  to  consider  large  displacements. 
For  example,  a  careful  snap-through  investigation  of  a  spherical  cap  demands  an  element 
allowing  for  varying  curvature.  Naturally  all  integrals  have  to  be  evaluated  numerically,  but 
this  presents  no  problem  with  a  set  of  suitable  pivotal  points. 

Note  that  for  the  adopted  idealization  the  definition  of  membrane  state  and  corresponding 
nodal  displacements  is  exactly  the  same  as  that  for  bending  in  the  TUBA  element.  Thus  we 


355 


AFFDL-TR-68-150 


effectively  prescribe  at  the  vertices  for  all  displacements  not  only  their  values  but  also 
their  first  and  second  derivatives.  In  addition  the  first  tangential  derivatives  normal  to  the 
edges  are  introduced  for  all  displacements  at  three  intermediate  points  on  the  edges. 

The  theory  is  most  easily  established  by  using  a  judicious  combination  of  vector  and 
matrix  calculus.  Scalar  and  vector  products  appear  in  this  connection  and  are  indicated  by 
the  standard  symbols;  e.g.  Q  -b  and  a  »b  respectively.  Moreover,  in  order  to  avoid  any 
misunderstandings  the  connotation  vector  is  restricted  in  the  present  note  to  true  vectors 
and  not  applied  to  column  matrices.  It  should  be  emphasized,  however,  that  occasionally 
matrix  operations  are  carried  out  with  actual  vectors. 

The  present  work  is  concerned  with  thin  shells  where  transverse  shear  deformations 
may  be  Ignored.  When  these  become  important  ASKA  uses  at  present  special  versions  of  the 
LUMINA  (Reference  3)  and  HERMES  (Reference  4)  sets.  Due  to  lack  of  space  the  derivation 
of  the  kinematically  consistent  load  and  mass  matrices,  as  well  as  of  the  geometrical 
stiffness  has  to  be  omitted.  These  are  discussed  in  Reference  19,  together  with  some  other 
interesting  features  of  the  theory. 


2.  DIFFERENTIAL  GEOMETRY  OF  SURFACES  IN  HOMOGENEOUS  TRIANGULAR 
COORDINATES 


a.  Definition  of  Geometry 


Let  us  assume  for  the  time  being  that  the  middle  surface  of  a  triangular  shell  element 

12  3 

is  defined  with  respect  to  a  global  cartesian  system  Ox  x  x  by  the  coordinates 


=  x'  ,  £  2>^3  *  '»  *2  =  x2(  ’^2  ’  £3  5 


3  3 

X  =  X 


(4a) 


which  are  functions  of  (i  =  1,  2,  3),  the  homogeneous  triangular  coordinates  in  the  basic 
triangle  1,  2,  3  (Figures  22  and  23).  The  functions  4a  give  a  unique  mapping  of  an  equilateral 
triangle  with  the  height  ‘one’  in  a  Euclidian  plane  onto  the  triangular  curved  middle  surface. 
The  position  vector  to  a  point  on  the  middle  surface  may  now  be  written  as 


x  =  {  x1  x2  x3  }  =  x  (■£,  ,  C2  .£3  ) 


(5) 
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Figure  22.  SHEBA  6  Shell  Element 


Figure  23.  Natural  Definition  of  Middle  Surface 
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b.  Metric  Relations  on  the  Surface 


The  incremental  vector  in  the  tangential  plane  is 


dx  =  x>(  d£t  +  x,2  dC2  h  x,3  dC3 


Figure  23,  where 


d  x  _  f  d  x1  dx2  dx3 

'i  =  d5i  '  1  acr  ac.  ac. 


The  incremental  arc  length  ds  is  obtained  from 
ds2  =  dx  •  dx  =  d  x2 

-Xj|  d£  (  +x,2  d£2  +  x  >3  d£3  +■  2X  I(  •  xi2  d£(  d£2—  2*,  2‘  *  ,3  d£2d£3+  2*  l3  *)(  d£3  d£i 

(8) 

Using  the  obvious  relation 

dC,  +  dC2  +dC3  =0  (9) 


one  establishes  easily  the  matrix  expression 


10  I  dC,  dC2 

i  i  O  dC2d?3 


0  I 


d£  d£ 

b3  *1 


Applying  Equation  10  in  Equation  8  we  find 


-ds  =(*,2-*,,)  d£,  d£2  +  (*,3 -xi2)2d^2dC3  +(».,  -x,3)2  d^dC, 

Ih  what  follows  we  denote  a  curve  on  the  middle  surface  as  a 

a  U  d£3 

/3  —  curve  if  £,  is  constant  ond  therefore  d£,  is  zero 


From  Equation  11  one  verifies  that  a  typical  incremental  arc  length  dsQ  on  an  a  -curve  is 
given  by 

dsa  )2d£, 


Since  dT^  =  0  it  follows  also  from  Equation  9  that 


dsa  --ix-z  -  <  =<«.,  )2  < 
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Note  that  d£^  and  d^  are  n°t  independent  in  Equation  13.  We  observe  that  the  arc  measures 
on  the  a,  (3 ,  y  -  curves  are  now  defined  by  the  moduli 


m 

a 


(14) 


c.  Natural  Derivatives  of  First  Order 


Following  our  usual  practice  we  introduce  the  operators 


a _ a  |  a  .  a  a  a  . 

dia  '  as,  ac2  •  a^‘~nz  <5C3  1 

As  in  Reference  12  we  also  apply  the  abbreviation 


(15) 


(16) 


where  fj  takes  the  values  a  ,  /3,  y 
Equations  14,  15,  and  16  that 


and  \f/  is  any  vector  or  scalar  function.  It  follows  from 


We  denote  as  natural  derivatives  the  total  derivatives  along  the  a-,  /3~,  y  -curves  with  respect 
to  the  corresponding  arc  lengths .  For  a  typical  /j,  -curve  we  have 

d  _  dgi  a  dCj  a  _  _j_  a  ^)8) 

ds^  ds^  aCj  ds^  aCj  a£^_ 


The  derivatives  dC.  /  ds  are  obtained  from  Equation  13 
k  H- 

fixed  in  accordance  with  the  conventions  in  a  plane  triangle. 


and  the  signs  of  the  roots  are 


The  natural  derivatives  of  the  position  vector  x  are  the  so-called  ‘natural’  unit  vectors 


e 


a 


d  x 

ds. 


*0  s 


d  x 
d  sy 


d  x 


X 


ds 


(19) 


'a  ~5(3  "y 

which  are  lying  in  the  tangential  plane  to  the  middle  surface;  see  Figure  24.  The  normal  unit 
vector  to  the  middle  surface  is  obtained  from  the  three  equivalent  expressions 

xe. 


ea  xe/3 


a 


x  e 


e/3 


x  e. 


JL 

x 


•y  x«a 


VXea 


(20) 


Following  Figure  24  the  moduli  in  Equation  20  are 


ea  x  e/3 


=  s  i  n 


X 


e/S 


x  • 


X 


:  s  in  a 


V  xea 


=  si  n  /3 


(21 ) 
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loco!  incremental  triangle 


a  -curve 


tongent iol  plane 


Figuxe  24.  System  of  Natural  Unit  Vectors 


where  a  ,  J3 ,  y  are  here  the  angles  in  a  local  infinitesimal  triangle  formed  by  the  a,  (3,  y 
curves.  The  distance  along  the  normal  n  is  denoted  by  z.  It  is  occasionally  convenient  to 
introduce  the  unit  tangential  vectors  fQ  f  t  f  ^  normal  to  eQ  respectively  and 

defined  by 


fQ  =  n  x  «Q  ;  f/3  =  n  x  •£  ',  ty-.  n  x  ey 


(22) 


(see  Figure  24).  ft  follows  from  Equation  15  that 

a  d 


and  hence 


dt/3  dSy 


=  0 


(23) 


x,  +  x,n  +  x,  ;  0 

a  (3  r 


(24) 


in 
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One  therefore  obtains  from  Equations  18  and  19  the  useful  relations 


ma  —  +■  my3  ~T —  +  mv  —r- —  =  0 


ds, 


£ 


y  ds 


r 


ma  ea  myS  *0  +  mx  ®7  =  0 


(25) 

(26) 


For  typographical  brevity  we  sometimes  denote  the  total  derivative  of  a  vector  ^  with 
respect  to  a  typical  natural  arc  length 

d 


ds^  r’""> 


When  applying  this  relation  to  the  a,  -curves,  m  takes  in  turn  the  values  a,  b,  c. 


d.  Natural  Derivatives  of  Second  Order 


[27) 


Using  Equation  18  one  finds 
^2 


V 


d  /  d  \  2 

ds^  (>  ds^T'V 


dm 


ds£ 


'/i.  d 


dv 


We  confirm  with  Equation  17  that 


dm 


I 


/.  a2x  \ 

v  ar.2  )  • 


ds/A 

Equation  28  yields  now 

j2 


m 


H-  ■  * 

i 


"ST  ‘V  •  ‘W 


ds. 


m. 


m  (  V  'Vm5  d£t 


(28) 


(29) 


(30) 


Application  of  Operator  30  on  x  and  use  of  an  evident  extension  of  the  Convention  27  to  second 
derivatives  leads  to 


d2  x 


ds 


,  =  x 

2  ’mm 


>  ft 

Equation  31  shows  immediately  that 


d  dp. 

ds..  =V’m 


I 


m. 


[  *' fi.fi.  (  V  V  ] 


H-  'H-H-  M 


de , 


ds. 


-  0 


which  follows  also  from  •^  =  I  .  Similarly  one  obtains  from  n  =  I  and  n  •  -  0 

d  n 


dn  den. 

n  •  — r =  0  and  n  - 


ds. 


ds , 


(31) 


(32) 


33) 
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In  accordance  with  the  differential  geometry  of  surfaces  the  radius  of  curvature  in  a  normal 


plane  n  -  is  given  by 


_J _  dn 

dv*C/A 


n  '  *>mm  "  mz^  n* 


e.  Variation  of  an  Element  of  Arc  Length  with  z 

Figure  25  shows  that  the  position  vector  x'  to  a  point  at  a  normal  distance  z  from  the 


middle  surface  is 


x  =  x  +  z  n 


Let  us  write  the  incremental  tangential  vector  dx  between  two  points  x,  and  x2  on  the 


middle  surface  as 

dV  "*2  “  Xi  =  v  dv  (36) 

Then  d  it  *  in  a  parallel  surface  at  a  distance  z  is 
r0 

it'  =  x'2  -  x'(  =  (x2  +zn2  )-(x|+zn|  )--dx^  +  z  dS/x  =(e^  +  z  <37) 

fL  jj. 

where  the  differential  form  arises  from  Taylor’s  expansion.  Neglecting  terms  of  higher 


order  one  obtains 


.2,2  dn  2 

(dx^)  =  {  ds^t)  =  (l-HZze^  •— — )(ds^) 


or  to  a  first  approximation 


f.  Element  of  Area 


d‘A  =  ‘ 1  +  ZV'~^  ldsM 


To  establish  a  theory  of  thin  shells  we  have  to  consider  here  only  an  element  of  area  dXI 
on  the  middle  surface.  On  a  a  -curve  we  have  e.g.  an  incremental  tangential  vector 

dxa  s  XM  d£i  +’x*2dC2  =<*,2  -  *„  >  dC2  =  x,a  d^2  (40] 

and  similarly  on  /3-and  y  -curves 


dX/3  =  X'/S  d^3  ;dxr  =  *•/  ^ 
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Figure  25.  Dependence  of  Arc  Length  on  z 


Now  d  vQ,  may  be  written  as 

dfl  =  |d»a  x^lHd-5  *"«yl  ■  |  dxyxdIa  I 

1  I x  ’a  *  x  ’/?  I  d^3 :  I x  ’/3  x«.y  K  <,  =  lx-y  xx-alde.  «z 

Forming  the  vector  product  of  %  with  Expression  24  we  confirm  that 

s  s  I  x,a  xx-/8  Is  I  x,$  x  x,y  I  =  lx.yxx.a  I 

Using  Relations  18,  19,  and  21  we  find 

S  =  m  m n  siny  =  m0  m  sina  =  m  m  sin B 
a  p  '  p  y  y  a  ^ 

3.  STRAIN-DISPLACEMENT  RELATIONS 
a.  Definition  of  Displacements 


(42) 


(43) 


(44) 


Let  us  assume  that  the  global  cartesian  displacements  of  the  middle  surface  are  given  by 
the  functions 


u'.s  u'(£ i  .£ 2  .£ 3  )  u2  =  u2  (C  , t  ,C  ) ;  u3  * u3  (C.  , C2  ,£,  ) 
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We  now  define  the  displacement  vector  of  a  point  in  the  middle  surface  as 

u  r  |  u1  uz  u3}  =  “  (£,  ,  C2  >^3  1  (46) 

As  mentioned  in  the  introduction  it  is  necessary,  in  order  to  develop  a  consistent  theory,  to 
express  the  geometry  and  displacements  in  an  identical  manner.  This  is  in  contrast  to  standard 
shell  theory  which  usually  operates  with  displacement  components  in  the  principal  directions 
of  curvature  and  in  the  normal  direction  to  the  surface. 

b.  Natural  Membrane  Strains 

Following  Figure  26  we  consider  a  typical  element  dx,,  of  a  -curve  between  points 

r’ 

x  (  and  x£.  Using  a  bar  to  denote  a  deformed  state  we 

Ti  =  X,  +  u,  ;  T2  =  x2  +U2  (47) 

Following  an  argument  similar  to  that  in  “Variation  of  an  Element  of  Arc  length  With  z”. 


Figure  26.  Total  Natural  Membrane  Strain 
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and  therefore 


dV  =  <  1  +  V  U*  m  )dV 


(49) 


The  total  natural  strain  in  the  middle  surace  and  the  /j.  -direction  is  now  by  definition 

du 


ds/i  “ds^ 


ds . 


>  ds,; '  V  'u>m 


(50) 


afj.  r~ 

In  the  last  expressions  the  shorthand  notation  du  /d  =  U,mof  Equation  21  has  been  applied. 
The  matrix  of  the  natural  membrane  strains  is  therefore 


N 


(51) 


Consider  next  an  element  dx,,  at  a  distance  z  from  the  middle  surface.  If  we  neglect  for  the 

r* 

time  being  the  change  of  the  normal  due  to  deformation  we  observe  that  the  total  elongation 
of  the  element  dx^  is  the  same  as  i 
is  modified,  however,  and  becomes 


of  the  element  dx',,  is  the  same  as  that  of  dx,, .  Since  the  arc  length  is  z  dependent  the  strain 
r-  r1 


ds  i 


I 


,  ^  _ 

"  H-0  d«^  ■  I  +  (z  /Rpfj) 


(52) 


where  the  last  relation  follows  from  Equation  39. 


c.  Rotation  of  Natural  Bi-peds  due  to  Deformation 


Let  us  consider  the  displaced  ‘unit’  infinitesimal  arc  elements  in  the  a-  and  /3  -directions. 
Following  Figure  27  we  have 


•+ 


d  u 


a  a  r  ds. 


—  *  du 

’0  “0  +  "dS 


(531 


The  unit  normal  vector  in  the  deformed  state  is  then 

_  »a  x'i8 


6  a  x  */3 


154) 


Applying  Equation  53  and  neglecting  higher  order  terms  in  the  derivatives  of  the  displacement 
vector  one  finds 


ea  x«/3  =  ea  xe£  +  xu’b 


u-a  xay 


(55) 
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or 


•a  *•[)  >' 


'a  x,/3 


.  Z 

-  sm  y 


)Z+2  CQ 

x*/3  '  x“'b  +% 

,  .  2 

"■••a  Xu'b  +  "'a 

t  1 

sin  y 

It  follows  that 


'ax«/9| 


sin- 


['  "  ’“a*  +  “.„x<*3  »] 


(56) 


(57) 


Substitution  of  Equation  57  in  Equation  54  and  application  of  Equation  55  yields,  after  some 
vector  transformations. 


m  =  n  —  n  = 


s  in  y 


[n  x(  aQ  xu.b  +um  x*(3  >]  x  n 


(58) 


where  m  is  the  change  of  the  unit  normal  vector.  Expression  58  may  be  interpreted  as  a 
rotation  of  n  through  the  angular  vector  ,  i.e. 


m 


=  x  n 


with 


a. 


sin 


X 


(59) 


(60) 


Using  Relation  25  applied  to  u  as  well  as  relation  26  it  may  be  readily  proved  that  the 
following  alternative  formulations  hold  for  Si  n 

fl"  --7^  [<n  '  W  =  FH5  [«"  •  V'r'"  •  VV>  1555  [<"  •  <61> 

To  obtain  the  angular  vector  of  the  bi-ped  n-e^one  has  to  add  to  £2n  a  rotational  vector  i2„n 
parallel  to  n  which  expresses  the  missing  rotation  of  about  the  normal  n  .  Clearly, 
following  Figure  27 


=  [  u,m  -<"*  V’]  "  =<  *'m  '  V  “ 


(62) 


The  total  rotation  of  the  bi-ped  n  -  e^is  thus 


&u  z  + 


(63) 
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d.  Natural  Bending  Strains 


Figure  28  shows  that  the  vector  d  x^  In  a  parallel  surface  at  a  distance  z  Is  transformed 


due  to  the  rotation  vector  SI..  Into 

H- 


dx'  -  d  x '  -I-  ii  x  d  x '  +  z  (^ds  XnJ 


ds/x  ^ 


(64) 


Since  the  second  and  third  terms  on  the  right  hand  side  are  small  quantities  n2  may  be 
replaced  by  n  .  We  next  find,  using  the  identity  dx^-  (  Sl^  *  dx^  )  =0 

( d^V  \  > 

(-^Tlxn)dVdV  ,65> 


>  )2s(dv'z 


,  .2 

> 


(d*„  )  :  (ds„  )  =  (  ds„  )  +  2z«  .(  - —  Xn]ds,’,  ds, 

h-  V  d! 


Hence  we  deduce 


d  V  =  dV  +  z  V  xn  )dv 


(66) 


Defining  in  the  usual  manner  the  measure  of  bending  strain  by 

X„ 


d  s^.  —  ds 


z  ds , 


(see  also  Reference  31),  we  obtain 


(67) 


(68) 


For  orthogonal  coordinates  this  expression  can  be  compared  with  Gol’denveizer’s  (Refer¬ 
ence  31)  corresponding  vector  expressions.  The  actual  bending  strain  at  a  distance  z  is 
therefore  given  by 

rl  c  ^ 

Z 

X, 


^  ds^  d 


h  ds^  l  +  (z/R 

We  form  the  three  components  Xq  ,  X^g,  X  as  the  column  matrix 

*N  1  X/9  Xy) 


(69) 


(70) 


e.  Formulae  for  the  Bending  Strains 


We  consider  here  the  typical  bending  measure  XQ  ,  which  is  split  into  two  parts 

d  id.  d  Shr 


Vxon  +x 


:  n  •  (  e  x  - —  \ 

aa  \  a  ds a  i 


an 

ds.„ 


(  n  x  eQ  ) 


(71) 
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Let  us  first  derive  dXln /dsa  .  From  the  first  formulation  of  Equation  61  we  find 

V/3  ]  - 

Writing  Expression  21  in  the  form  sin  y  =  n  •(  e^x  e^i  we  deduce,  using  the  first  of  Equation  33 


J  sinX  /  d«q 

dsa  dsa 


xeo  1-  ea  x 


Substitution  of  Equation  73  into  72  yields  a  lengthy  expression.  Carrying  out  however  the 
vector  operations  indicated  in  Equation  71  one  finally  obtains  using  the  Convention  27  for 


n,  u  ,  e. 


y  =  -  n,  u  -  n  u 
Aan  a  ’a  'aa 


+  ^"-K,o''a  |u-b  ‘‘Vo  V'o  J  (74> 

In  order  to  symmetrize  Equation  74  we  carry  out  corresponding  operations  on  the  third  formula 
in  Equation  61  and  form  the  mean  value  of  the  two  equivalent  expressions  to  find 

Xan  *a  *a  *00 

-tKo  1751 

For  the  second  term  XQ  aone  readily  derives 

W  ‘ "  XV  >]  [-^  <"  *V  ']=<».,•»«  IK„  la  1  176) 

Combining  Xan  and  XQa  one  finds  the  relatively  simple  formula 

Xa  =  *a,a  >'V  “'0  1  -  "  ■“•oo 

"tK,o  V)]‘"'u'0  l  +  i'1Vo’,a,["'(iy  U'b-I^5  v)]  1771 

Corresponding  egressions  for  X ^  and  X  are  easily  established.  For  lack  of  space  we 
have  to  omit  the  proof,  that  rigid-body  motions  indeed  give  neither  membrane  nor  bending 
strains. 

4.  INTERPOLATION  OF  THE  MIDDLE  SURFACE 
a.  Definition  of  Geometry  at  the  Vertices 

It  can  be  assumed  that  the  lines  of  principal  curvature  at  the  nodal  points  are  known  and 
that  their  directions  are  given  by  the  unit  vectors 


A  ’  dA 


8  '  dB 
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where  the  derivatives  are  taken  with  respect  to  the  arc  lengths  on  the  principal  curves,  in 

the  case  of  cylindrical  shells  or  shells  of  revolution  it  is  no  problem  to  find  the  e  “  e_ 

A  B 

system.  For  a  general  surface  we  may  apply  the  same  mapping  procedure  as  for  the  LUMINA 
element  (Reference  3).  To  this  purpose  we  only  require  the  cartesian  coordinates  of  a  set  of 
nodal  or  dummy  points  to  establish  a  representation  of  the  shell’s  middle  surface  by  the 
Lagrangian  interpolation  polynomials  v  (  £  £2)  .  Following  the  solution  of  a  simple  (2  x  2) 

eigenvalue  problem  it  is  easy  to  transform  the  derivatives  of  x  (£  '  £  2)  with  respect  to  £  1 
and  £  to  the  so-called  principal  set 


XX  XX  X  X 

»A  ’B  ’A  A  ’AB  ’BB 


(79) 


For  the  transformation  of  the  second  order  derivatives  we  have  to  apply  a  technique  corre¬ 
sponding  to  that  of  Equations  92  and  95  for  specific  details  see  Reference  19. 


b.  Definition  of  Geometry  within  an  Element 


To  obtain  the  necessary  data  defining  the  geometry  within  an  element  we  have  to  transform 
the  nodal  data  79  into  the  natural  directions  of  the  triangular  shell  element.  Following  Figure  29 
we  first  introduce  the  unit  vectors  of  the  basic  triangle  i ,  j  ,  k 


J, 


*j  “  *j 


Xi,  —X 


X 


J, 


X  ;  -  X  , 


80) 


A  '  -> 

We  next  determine  at  a  typical  vertex  i  the  natural  unit  vectors  e.  and  e  .  using  the 

Ai  v\ 

conditions  that  e^.  (  e  . )  lie  in  the  plane  nr  j^(  nr  J  )  and  also  in  the  tangential  plane.  We 
easily  find 

i\  >  " 


ex  = 


[i-  (njx)2  ] 


Vz 


(81) 


where  the  subscript  i  is  omitted.  The  actual  interpolation  procedure  is  carried  out  in 
Equation  98  using  the  TUBA  6  functions.  To  be  consistent  with  the  TUBA  theory  it  is  necessary 


to  fix  the  measures  at  the  vertices  i,  j,  k  using  equations  of  the  type 

d  x  _  j_  dx  _  _j _  _  i 

dtx  "  /X  >X  ’  \  -  i  \  >  X 


(82) 


where  tx  is  a  distance  measured  parallel  to  jx  in  the  basic  triangle;  see  Figure  29.  On  the 
other  hand  dx  t  dtx  has  to  be  identical  with  ex  •  Therefore  we  have  to  set 


\~ - ; — 8.  and  correspondingly 

A(e,*K)X 


(•»  *  K  ) 


(83) 
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which  yields  for  the  associated  measures  at  the  vertex  i 

mx  “A  ’  '■  "V  =  •  K  ' 

To  obtain  the  last  at  the  vertex  i  we  apply  Equation  26  to  find 


(84) 


V  =-^rr(m\eA  +"V  «i/  ) 


and  since  I  e/f  I  =  I 


=  (ml  +2m,  m  (*  •  e  )+mz  )/z 
r  K  A.  "  X  v  v 

Using  the  unit  vectors  eA  and  e  g  of  Equations  78  we  find  the  evident  identity 

eXs  { *X  '  ®A  *  *’A  +  (  W“’B 


(85) 


(86) 


(87) 


and  a  corresponding  one  for  .  Equation  83  yield  now  the  transformation  of  first  order 
derivatives 


*u  ■  "’X<*X  ,A,X'A  mx*  *  X  '  VS 

Vs  W  V*-A  +  "V  ‘V  V  “'B 

which  is  used  to  define  the  operator  transformation 


(88) 


where 


IJ-  d 


dA  dB 


F. 


m .  m 
1  X  v  J 


•X'* A  V*B 

*v  “b  J 


e  .  e. 

v  A 


(89) 


(90) 


To  extend  the  Transformation  90,  which  is  exact,  to  second  order  derivatives  we  assume  it  to 
represent  a  transformation  from  an  oblique  system  £  -  £  to  an  orthogonal  one  A  -  B  in  the 
tangential  plane,  ignoring  any  change  of  the  unit  vectors  eA  and  eB.  This  leads  to 


where 


V 


* 


a2  a2  i 

5  r  a2  a2  a2 

1  «x 

• 

1 — » 

w  * 

& 

X_A 

'  <•  dA^  dA  dB  dBz  J 

»x  -.A  I* 

2(e^eA)(e  eB) 

"W  J 

(*X  vv  -A1  <*x 

VW+'YVW 

(91) 


.2  n 


<v V 


2(V*A  )(V*b  > 


<*x  *b> 
YV(VB 

<•*  V 


(92) 
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However,  it  follows  from  Equation  23  that 


a2 

a2  , 

„  a2  __  .  a2 

(93) 

aC4  '  «x2 

Thus,  one  may  replace  Transformation  91  by 

r  a2  a2  a2  •»  _c 

aCx  K  < 

,  d‘  a2  a2 

\ 

(94) 

1  ^A  Ku 

where  '  r 

1 

< 

%  i 

1 

1  aA2  a*  aB  aBz 

2<*x  •a>(,a  ,b> 

J 

(eA  ‘eB  )E 

si  •  K  v  mi  J 

*  ®A  )Z 

zi  V  •,a)(V'*21 

(®/i.  '®B)2 

(95) 

2<*„  •b> 

In  establishing  S  (  Relations  85  and  86  have  been  used.  Applying  the  Transformations  89  and 
95  we  can  derive  at  each  vertex  from  the  principal  set  79  the  natural  set 


X  X’\  %  X’A\  %'vv 

c.  Interpolation  of  the  Middle  Surface 


(96) 


In  order  to  use  the  TUBA  6  interpolation  scheme  we  require  the  (21  x  3)  natural  geometry 
matrix 


IN 


=  {xi  *2  x3  xi,a  \r  \/3  X2,a  x3,/  X3,/Q  xi,aa  *1,77  X2,^S  X2,aa 


IN 


X 

3 

Xi,a  x 

»,j80  *ij8/3 

2 

x  3 

IN 

IN 

3,77  *,J3/3  *1  ,(3(3  2,77  3-aa  t’V*  5>I/I  *>vz 


} 


(97) 


The  elements  in  Equation  98  are  obtained  by  transposition  of  the  column  matrices  97  with  the 
exception  of  x  _  ,  x  _  and  xe  , . ,  which  are  derived  by  the  interpolation  procedure  of 

Equation  39  in  Reference  12.  Note  that  the  eccentricity  parameters  /i  refer  to  the  basic 

points  4,  5,  6.  The  alternative  formulation  of  xJN  in  Equation  97  represents  it  as  a  matrix 
of  three  columns  corresponding  to  each  of  the  three  cartesian  coordinates.  In  accordance  with 
Equations  6  and  24  of  Reference  12,  we  find  the  current  position  vector  x  from  the  trans¬ 
formation 

xf  =  x  =  u)  A  x  IN  (98) 


The  (1  x  21)  matrix  <0  which  depends  on  £2,  and  the  (21  x  21)  matrix  A  are  given 
in  Equation  31  and  Figure  5  of  Reference  12.  The  transformation  98  may  be  interpreted  as 
defining  three  distinct  surfaces  £2»  C3)  over  the  basic  triangle.  The  geometrical 
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data  associated  with  adjoining  triangular  shell  elements  meeting  at  a  nodal  point  are,  however, 
deduced  from  a  common  smooth  surface.  It  may  hence  be  argued  that  the  functions  x  P  and 
their  first  derivatives  are  continuous  along  the  edges  as  in  the  case  of  the  deflexions  w  of 
TUBA  6.  5  is  of  no  consequence  that  the  basic  triangles  are  not  in  a  plane  but  form  a  polyhe¬ 
dron  surface.  Consequently  we  have  interpolated  the  actual  middle  surface  of  the  shell  in  such 
a  manner  that  no  folds  are  generated.  The  second  derivatives  are  of  course,  in  general, 
discontinuous  across  the  edges  with  the  exception  of  the  vertices.  We  may  also  define  a 
(63  x  1)  matrix 


IN 


■{ 


!,0£  I,/ 


X 


3,aa 


4,  J/3  5,  V\  6,1/2 


(99) 


which  yields 

["  A]  X1N  *  D,IN 


(100) 


where  the  operator  3  is  defined  in  Reference  5. 


5.  INTERPOLATION  OF  DISPLACEMENTS 


a.  Displacement  Matrix  at  the  Nodal  points 

At  a  typical  vertex  i  of  a  triangular  shell  element  i,  j,  k  we  introduce  the  (6  x  3)  dis¬ 
placement  matrix 

fi,  '■  {  “  “'A  “'B  **’  AA  ■•AB  %B  },  1101 

Applying  the  Transformations  89  and  94  we  deduce  the  natural  displacement  matrix  at  the 
vertex  i 


with 


002) 


003) 


At  the  intermediate  nodal  points,  p  ,  q  ,  r  (see  Figure  29)  we  have  to  define  a  (1  x  3)  dis¬ 
placement  matrix  in  accordance  with  the  TUBA  procedure.  For  instance,  at  the  point 
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The  differentiation  with  respect  to  u.  is  fixed  by  operators  of  the  type 


dv 


=  2 


-(af T  +  t-)-^( 


-+ 


H, 


■) 


105) 


■3  -  >3  '  I 

see  also  Equations  16  and  17  of  Reference  12.  No  further  transformations  will  be  applied  to 
these  freedoms. 


b.  The  Complete  Displacement  Matrix  of  the  Element 

We  next  introduce  the  complete  (21  x  3)  p1  matrix  for  a  shell  element  in  the  form 

P  z  { P  P  P  P  P  P\  (l06) 

X  '•  I  2  3  4  5  6 

where  p  ^  pz ,  are  given  in  Equations  102  and  p* ,  p$,  in  Equation  105.  The  complete 

natural  displacement  matrix  is  given  by 


^  ftl  ft 2  ft 3  ia  r4  l 


NI 


4ii  -  -  4fl  r  i 

*  T'=  V‘} 


Extending  the  Transformation  103  we  find 


P  =XS  p 

NT  s  1 


with 


(107; 


:  108) 


XS  =B.  XB 

where  B6  is  given  in  Equation  31  of  Reference  12  and 


Also 


4  Sit 


-i 


(21  X  21) 


*  *k  h 


:  io9) 


110) 


(!!!) 
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Alternatively  we  can  define  the  two  (18  x  1)  matrices  of  displacements  at  a  nodal  point 


*  -{ “ 


u  u  u 

’A  ’B 


p  ={  u  u,x 
Ni  1  A 

and  the  (3  x  1)  matrix  at  nodal  point 


'v 


AA 

U,AB 

U,BB 

XX 

U, 

% 


p,  n  k 


These  form  the  alternative  (63  x  1)  complete  displacement  matrices 

P  P  P  P  P  P 

I  ^  r|  '2  3  4  5  6  J 


and 


^NI  "  i  i  ^N2  3  ia  P 

The  Transformation  109  becomes  then 


4ft  _  4ft  4ft 


^(3  P* 


^Nl"  31  XS  "  Vi 


( 112) 

(113) 

(114) 

(115) 

(116) 

(117) 


The  conditioning  of  the  stiffness  matrix  may  be  improved  considerably  by  introducing  the 
displacement  matrix 


p  ■-  f  C 
rl  I  le 


c,  c 

ze  3e 


1  P 

9  J  ~ 


-  C,  P 


where  a  typical 


and 


=[  c. 

ie  1  i 


c.  c.  c.  c.  c 


,  J 


c. 


r 

€ 


■], 


(118) 


1119) 


(120) 


This  transformation  relates  the  global  cartesian  displacement  u  - {  u1  u2u3|  to  the  local  system 
principal  curvatures.  For  brevity  we  restrict  ourselves  in  the  above  derivation  to 
shallow  shells.  The  extension  to  more  general  surfaces  is  straightforward,  see  Reference  19. 

c.  Interpolation  Scheme  of  Displacements 

As  in  the  case  of  the  interpolation  of  the  middle  surface  the  current  displacement 
vector  u  is  defined  by  the  transformations 


ti  =  Ui  A  p 


NI 


or 


U  ;  Dp 


NI 


121) 
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where  D  is  given  in  Equation  101.  Note  that  the  (21  x  3)  matrices  of  a  pure  translation 
UQ  or  a  rigid  rotation  through  a  vector  <f>  with 

uo  s  {  uo  uo  uo  }  >  £  =  {  4>'  4>Z  <t>*  }  {l22: 

have  to  be  defined  as 

=(“o  "o  “o'  °,e,3  }  ('23: 

or 


NI<jb 


•[*' 


2  3 
*NT 


.3  2 

-(f)  X 


NI 


,3  I 

+  ‘ni 


, }  3 

-(f)  X 


NI 


,i  2 


NI 


,2  I 

-<p  X 


NI 


(124) 


respectively,  where  x ^  are  the  (21  x  1)  matrices  which  form  the  columns  of  in 

Equation  98.  Using  the  Definitions  (123)  and  (124)  the  TUBA  6  interpolation  scheme  guarantees 
that  at  any  point  within  the  element  the  relations  U  =  UQ ,  U  =  x  x  hold  and  that  there  are 
zero  strains. 


6.  STIFFNESS  MATRIX 


a.  Virtual  Work 

If  an  element  of  volume  dz  dfi,  is  subjected  to  a  strain  the  virtual  work  due  to  a 

stress  E  K  €  is 
pv  v 

d3  W  dz  dft  (125) 

where  E  is  standing  for  Young’s  modulus  and  K^JLV  is  one  of  the  natural  material  stiffness 
coefficients  which  form  the  matrix  (Reference  1).  To  evaluate  Equation  125  we  can  apply 
the  formulae  of  the  standard  shell  theory  (References  31  and  32)  to  the  direct  component 
strains.  We  obtain  after  integration  over  the  thickness  t 

d2  W  =  EK  ( t  €  €  +■  — —  Y  dil  (l26) 

fMV  /AO  VO  12  */A 

The  variation  of  the  elements  of  arc  length  and  area  with  z  can  be  neglected  in  the  case  of  thin 
shells  with  t/R  «1,  R  being  a  typical  measure  of  the  radii  of  curvature.  It  is  worth  mentioning 
that  the  natural  stiffness  concept  leads  also  to  a  variation  of  KN  with  z,  due  to  a  change  of  the 
enclosed  angle  between  and  ,  but  this  influence  is  once  more  of  the  order  of  magnitude 
t/R  (Reference  19). 
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b.  Membrane  Stiffness 


Using  the  interpolation  scheme  of  Equation  121  we  are  able  to  represent  the  Relation  52 
for  the  natural  membrane  strain  matrix  as  the  linear  transformation 


*No  =  m  °M  ^NI 


Here 


and 


m 


=  f  ma  myS  "V  1 


aM  =  m  { ea  D-a  */3  °>/3  Dv}  (3X63) 


where 


V,TI  [<v  A] 

The  membrane  stiffness  corresponding  to  the  displacement  matrix  p  is  now 

'NI 

V  E/n.  '  Vm"11  «63  * 


Here 


=  m"2  K..  m~z 


(127) 


(128) 


(129) 


(130) 


(3X3) 


(131) 


(132) 


Naturally,  the  Integral  131  and  all  following  corresponding  expressions  have  to  be  evaluated 
numerically. 

c.  Bending  Stiffness 


Introducing  the  Interpolation  Scheme  121  Into  Equation  77  we  find  the  linear  transformation 


.-2 


where 


V  m  o  p 

*N  B  PNI 


{  ° 8a  aB/3  aByJ 


(133) 


(  3  X  63  ) 


(134) 
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Restricting  our  considerations  to  the  typical  matrix  o  and  omitting  the  subscript  Ba  on 

60( 

the  righthand  side,  we  may  write 


a_  =  a  +a  +  a0  +  a 
Ba  aa  a  (3 


r 


(  I  X  63  ) 


(135) 


We  easily  establish 


where 


Also 


aaa  z  n  a  a 

1  bl  D 


1  b*  D. 


a  '  a  ’a 
■  -  '  »V  0 ,y 


(3-  ^  •  /  my 


ba  *  *  x,aa  n  ^  ea  2  [x>aa  (sin y  /3  sin  / 3  V)  ]  n 

'  V  =-r^g-U-aa-'a 


(136) 


(137) 


*’aa  ‘  D’aa  XNI 


(138) 


The  bending  stiffness  corresponding  to  the  displacement  matrix  p  is  therefore 

k  =  ~4~  f  *3  a*  XM  a  dfl  (63  X  63)  (139) 

B  12  B  N  B 

d.  Complete  Stiffness  Matrix 

The  complete  stiffness  matrix  based  on  the  displacement  matrix  p  (Equation  119)  is 
finally 

k  =  CI  XS  [kM  +  kB]XS  Ce  (63  X  63)  (140) 
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SECTION  IV 

ELASTOPLASTIC  ANALYSIS 

1.  GENERAL  CONSIDERATIONS 

The  matrix  displacement  method  covers  not  only  linear,  elastic  (Hookean)  behaviour,  but 
also  in  conjunction  with  step-wise  procedures,  nonlinear  conditions  which  occur,  for  example, 
in  nonelastic  and  large  displacement  or  strain  problems.  This  chapter  presents  an  extension 
of  the  theory,  first  given  in  Reference  36,  for  handling  nonlinear  stress-strain  relationships 
due  to  plastic  deformations  and  work  hardening.  In  order  to  limit  the  scope  of  this  paper,  we 
do  not  consider  the  effects  of  large  displacements  and  strains,  which  often  occur  in  such 
problems,  although  for  these  purposes  the  concept  of  geometrical  stiffness  has  been  developed 
(References  1,  6  and  Section  V). 

In  order  to  describe  clearly  the  principles  of  our  method,  we  restrict  ourselves  to  the 
von  Mises  yield  criterion  for  isotropic  materials  and  the  Prandtl-Reuss  equations,  which 
together  lead  to  a  mathematically  and  physically  consistent  theory  of  yield  phenomena.  The 
incremental  stress-strain  relations  corresponding  to  stepwise  small  changes  in  the  loading 
are  written  in  a  matrix  formulation  suitable  for  computer  use. 

As  pointed  out  in  previous  publications  (Reference  36)  it  is  advantageous  to  apply  the 
so-called  method  of  initial  loads  (Reference  1)  for  the  representation  of  plastic  strains.  This 
technique  allows  the  use  of  a  constant  stiffness  matrix  for  constant  temperature  and  small 
displacements.  The  calculation  of  initial  strains  and  the  elastic  analysis  with  initial  strains 
are  described  in  a  general  manner,  which  can  be  applied  to  any  element. 

The  determination  of  the  plastic  initial  strains  is  somewhat  complicated  because  the  in¬ 
crement  of  plastic  strain  is  dependent  on  the  stress  increment,  which  again  depends  on  the 
current  value  of  plastic  strain.  These  difficulties  can,  however,  be  circumvented  by  the  method 
of  this  report.  The  application  of  the  procedure  is  illustrated  on  the  analysis  of  a  cracked 
plate. 
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2.  INCREMENTAL  STRESS-STRAIN  RELATIONSHIPS 


a.  Stress  and  Strain  Vectors 

We  consider  only  a  three-dimensional  continuum,  where  strains  and  stresses  are 
described  using  the  <6  x  1)  column  matrices  of  cartesian  components 


t  £ 

{  €  €  €  -L-  €  4= 

l  xx  yy  zz  J2. 

€  ~~= 
yz  fz 

(142) 

cr  - 

\  cr  cr  a  -Jz  cr 

1  xx  yy  zz  v  xy 

y/2a 

yz 

(143) 

These  definitions  guarantee  the  duality  of  static  and  kinematic  transformations,  as  shown 
in  Reference  1.  Lack  of  space  precludes  giving  the  most  elegant  natural  formulation,  but  we 
suggest  for  this  purpose  a  perusal  of  Reference  23.  The  nonelastic  strains,  which  are  to  be 
understood  here  exclusively  as  plastic  strains,  and  the  total  strains  are  denoted  by  the 
vectors  17  and  y  respectively.  Thus  we  have 

y  =  €  +  17  ( 144) 

Clearly  17  and  y  must  be  written  analogously  to  c  of  Equation  142.  For  linearly  elastic 
material 

<r  =  E  «  (145) 


where  E  is  the  (6  x  6)  material  stiffness  matrix  (Reference  1),  which  in  the  isotropic  case  is 

E  =2S[16+  -Pfe  E3,3]  "46’ 

with  the  (6  x  6)  unit  matrix  I*  and  E,  =#__•*,  where  e,  ,is  the  (6  x  1)  column  matrix 

5  5,5  5,5 

e3  3  -  jl  l  IOOO| .  G  and  v  denote  the  shear  modulus  and  Poisson’s  ratio  respectively. 

The  vectors  a  and  y  can  be  divided  into  so-called  hydrostatic  and  deviatoric  parts 


<r  =  <r  4  tr 
H  +  D 


r  =  ru  +  yr 


where 


"■h  =ilo'.«+3-„  +  <r22|,3,3 '  r»  +  +rz2i*3,3 


(147) 


(148) 


are  the  hydrostatic  stress  and  the  delation  respectively.  To  fulfill  Equation  147  we  set  for  the 
deviatoric  vectors 


°-n  =  *  X  n  =  X 


(149) 
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where 


For  any  v  it  ia  easy  to  prove  that 

E^d  =  2G/i-D 

Note  also  the  useful  relation 

^D^D=^D 

for  the  symmetric  (6  x  6)  matrix  fj.  0 

b.  Yield  Criterion  and  Plastic  Strain  Increments 

Using  the  Definition  150  of  the  deviator ic  stresses  <rQ  and  Equation  152  the  von  Mises 
criterion,  which  is  valid  for  isotropic  materials,  takes  the  simple  form 

:  -  ff'  ff  :  -  (T  (153) 

*  2  D  0  2  Hi 

where  a  is  the  so-called  equivalent  stress.  Condition  153  for  the  stresses  <x  says  that  plastic 
deformations  occur  when  cr  reaches  a  critical  value,  which  is  initially  equal  to  the  elastic 
limit  o-^  .  Subsequent  yielding  takes  place  if  a  becomes  greater  than  the  maximum  value 
ct"  previously  reached  in  the  point  considered.  Corresponding  to  cr  we  now  introduce 
the  equivalent  strain 

yZ  =  y^  y^  =-| ~  y*  ( 154) 

For  a  uniaxial  stress  state  one  finds 

^  _  g.  .  y  -  y  .  rj  =  rj  =  y  -  cr/ E  (  E  =  Young's  modulus)  (155) 

One  sees  therefore  that  the  assumption  of  plastic  isotropy  allows  the  extraction  of  all  in¬ 
formation  about  the  yield  condition  in  a  three-dimensional  stress  state  from  a  simple 
uniaxial  tension  test. 

The  Increments  of  plastic  strain  7)^  of  a  typical  loading  step  are  calculated  using  the 
Pandtl-Reuss  equations 


(151) 


(152) 
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This  establishes  a  consistent  theory  in  connexion  with  the  von  Mises  Criterion  153.  Applying 
Equations  153  and  154  to  compare  the  three-dimensional  with  the  uniaxial  stress  state  we 
find  for  the  so-called  plasticity  factor  J3^. 


B  -11  °A 


(157) 


where  £  is  the  slope  of  the  cr-77  diagram  of  a  uniaxial  test,  as  shown  in  Figure  30.  All  un¬ 
loadings,  of  course,  must  be  calculated  elastically.  This  is  done  automatically  by  setting 
-1-  -  °i  ,  if  cr-crmaxor  .  The  equivalent  stress  increments  can  be  found  either  directly 
from  ex  =  or  by  differentiating  Equation  153.  This  leads  to 


°A 


cr1  a  <r. 
'd  A 


(158) 


where  cr ^  is  the  incremental  stress  vector.  Expression  158  is  especially  suitable  for  the 
following  theory  because  it  enforces  a  complete  linearization  of  the  incremental  elastoplastic 
stress-strain  relations. 


c.  Initial  Loads 


It  is  advantageous  to  represent  the  effect  of  the  plastic  strains  by  a  corresponding  set  of 

initial  loads.  These  simulate,  in  effect,  an  elastic  suppression  of  the  initial  strains  within 

each  element.  The  total  strain  can  be  obtained  for  any  type  of  finite  element  from  the  vector  p 
of  element  freedoms  by  a  transformation 


r  =  ~  p 


(159) 


At  the  same  time  the  initial  strains  of  any  point  within  an  element  may  be  determined  from 
a  characteristic  set  of  initial  strains  7)c  using  an  interpolation  transformation 

V  =  to  Vc  (160) 

where  77  c  is  a  super  column  matrix  containing  the  vectors  77  of  a  certain  number  of  suitable 
points.  The  virtual  work  principle  yields  now  for  the  initial  load  vector  j  associated  with  p 


where 


J  = 


tOu  E  77  d  V  = 


(161) 


h 


/  a  E  b  dv 

V 


(162) 
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Figure  30.  Typical  ct-tj  diagram  (above) 


3n  an  analogous  fashion,  a  set  of  total  strains  describing  uniquely  the  deformation  of  an  element 
may  be  set  up  in  the  form 


7c  P  (l63) 

where  /tu,  is  a  column  supermatrix  built  15)  from  the  matrices  xl  of  the  selected  characteristic 
points.  For  simplicity  in  what  follows  we  omit  the  subscript  c  and  do  not  introduce  ad¬ 
ditional  subscripts  or  new  symbols  when  expanding  matrices  defined  elementwise  to  complete 
structural  matrices. 


383 


AFFDL-TR-68-150 


The  condition  of  equilibrium  for  a  structure  subjected  to  external  loads  R  and  initial 
loads  J  -applied  to  the  elements  -  is 

R  =  Kr  -  o'j  ( 164) 


where  r  is  the  displacement  matrix  corresponding  to  R  ,  K  the  elastic  stiffness  matrix, 
and  a  a  Boolean  matrix  relating  the  degrees  of  freedom  of  the  elements  to  those  of  the 
structure  (Reference  1).  We  deduce  from  Equation  164. 


nK  '(R-fl'j  (165) 

The  total  displacement  vector  p j  of  the  elements  is  found  from  the  Boolean  Transformation 

/t>T  =  ar  066) 


which  yields  for  the  total  strains 


y  ~4V  pt 


067) 


Therefore  the  elastic  strains  are 

068) 

(169) 

d.  Incremental  Stress-Strain  Relationships  for  a  Complete  Structure 


€  =<0 'PT  -  7) 

and  the  stresses 

<r  =  E Ia/>t  -7)) 


Condensing  Equations  162,  166  to  169  for  a  structure  subjected  to  external  forces 
and  plastic  initial  strains  ij  ,  one  obtains  the  stress  increments 

<r^  =  E  A  a  K  '  +  E  (a  a  K~'  af  h  -  I  )  if  (170) 

where  I  is  a  unit  matrix.  Using  Equations  156  to  158  it  follows  that 


where,  in  agreement  with  Equation  170 


HR  + 

A 


H  =  2G  C  'SS^OK 


-! 


G  =2g£  SS  T 


with 


«■ 


T  =  A.Q  K-1  of  h  -  I 


(171  ) 


(172) 


( 173) 
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The  material  stiffness  E  has  been  eliminated  using  Equation  151.  All  slopes  £  of  the  cr-  7 7 
diagram  are  collected  in  the  diagonal  matrix  £  .  Note  that  by  definition  £  1  has  nonzero 
elements  only  at  such  points  where  further  plastic  deformation  occurs.  The  matrix  S  is  a 
diagonal  supermatrix  whose  elements  are  the  vectors 

S  =  —  —  M  n  O'  =  4  ~=~  n  (l74) 

2  <y  '  D  2  cr  D 

The  operation  SS*  represents  therefore  a  projection  of  the  elastic  strain  vector  €  onto  the 

normalized  vector  4-  <r .  We  also  observe  that  the  transformation  matrix  T  yields  the  elastic 
cr  D 

strains  <  due  to  any  initial  strains  V)  via  the  relation 

€  =  T  7)  (175) 

In  Equation  171  HR^  stands  for  the  increments  of  plastic  strain  calculated  from  stress  in¬ 
crements  entirely  elastic  and  arising  from  the  load  R^  .  This  contribution  to  7}^  can  thus  be 
easily  calculated  without  first  forming  H  .  The  second  part  shows  that  plastic  initial 

strain  increments  17  A  themselves  produce  further  plastic  deformation.  Therefore  it  is  nec¬ 
essary  to  know  G  explicitly  In  order  to  resolve  the  system  of  Equation  171  for  7 j 

3.  ITERATIVE  METHODS  OF  ELASTOPLASTIC  ANALYSIS 

a.  The  Basic  Problem  of  Iterative  Techniques 

The  methods  of  elastoplastic  analysis  described  in  Reference  36  use  estimated  starting 
values  7Ja_  for  17  A  on  the  right  hand  side  of  Equation  171.  This  procedure  enables  us  to  com¬ 
pute  improved  plastic  strain  increments  7) A+  without  knowing  G  explicitly.  To  examine  this 
point  further  we  rewrite  Equation  171  in  the  form 

•n  =  HRa  +  G  Tj .  (176) 

'A  +  A  'A- 

There  are  two  possible  methods  of  evaluating  Equation  176.  These  may  be  called  the  direct 
incremental  method  piM)  and  the  normal  iterative  method  (NIM).  DIM  uses  the  improved 
plastic  strains  17^  of  a  typical  loading  step  i  as  starting  values  7)^* 1 1  of  the  following  step. 
On  the  other  hand  NIM  iterates  within  each  loading  step  using  as  intial  17  the  results  of 
the  preceding  iteration  until  convergence  is  found.  This  can  be  understood  as  an  iterative 
solution  of  Equation  171,  which  can  converge  only  if  all  diagonal  elements  Gjj  of  G  are  less 
than  ‘one’  as  proved  in  standard  matrixtheory.  However,  this  condition  is  not  fulfilled  in  most 
practical  problems.  For  duct  He  metals  £  decreases,  in  general,  exponentially  with  increasing 
equivalent  stress.  Following  the  second  of  Equation  172  l/£  is  directly  a  multiplier  of  the 
diagonal  elements  of  G  .  Therefore,  with  increasing  load,  a  point  will  be  reached  at  which 


385 


A  F  FDL-  TR- 68-150 


any  one  of  the  G  -  j  may  become  greater  than ‘‘one”.  Especially  in  three-dimensional 
problems  it  is  even  possible  for  T  in  Equation  173  to  have  such  large  diagonal  elements  that 
the  iteration  does  not  converge  for  any  reasonable  value  of  £  .  It  is  now  obvious  that  the 
matrix  G  has  to  be  considered  in  order  to  develop  a  successful  method  of  solution  for  plastic 
problems. 


b.  Improved  Iterative  Method 

In  contrast  to  NIM  we  denote  the  procedure  which  uses  G  as  an  improved  iterative 
method  (VIM).  The  calculation  of  G  from  Relations  172  and  173  involves  extensive  numerical 
computations,  especially  for  the  complete  inversion  of  K  .  The  dimensions  of  G  ,  however, 
can  be  readily  reduced  by  a  factor  of  six  for  the  three-dimensional  case  by  introducing  the 
scalar  plasticity  factors  instead  of  the  vectors  7)^  in  Equation  171.  This  leads  to 


=  2 G  St  T  S 


(177) 


Obviously  a  numerical  evaluation  of  Equation  176  is  only  economical  if  an  overwhelming  part 
of  the  structure  is  plastically  deformed.  In  most  real  problems,  however,  yielding  occurs 
in  relatively  small  regions  near  notches  or  cracks.  We  require  then  only  a  certain  part  of  G  , 
obtained  from 

<5  =  C  G  Cf  (178) 


where  c  is  a  Boolean  (Nxn)  matrix  which  selects  the  N  characteristic  points  where  further 
yielding  takes  place,  out  of  a  total  of,  n  characteristic  points.  Rewriting  Equation  171  in  the 
functional  form 


’’a  *  *‘V 


it  is  easy  to  see  that  G  is  the  gradient  matrix  of  <f> ,  i.e. 


(179) 


G  = 


(180) 


Differentiation  with  respect  to  matrices  is  defined  in  Reference  1.  Since  it  is  not  difficult  to 
determine  values  of  for  any  ,  G  can  in  practice  be  replaced  approximately  by  difference 
quotients  computed  from  modified  vectors  7)^.  Let  us  assume  that  we  know  an  approximate 
value  V  /\q  with  N  plasticity  factors  P^ 0 .  Then  <p  is  calculated  not  only  for  tj  alone 
but  also  additionally  for  N  modified  vectors  Tjf  in  each  of  which  one  P^  is  slightly  changed. 
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From  the  differences  *  <f>  ( 77^- )  one  may  derive  gradients,  sufficiently  accurate  to 

ensure  convergence  of  the  iteration  described  below.  This  method  consumes  little  more 
computer  time  as  it  simply  requires  considering  N  +  1  loading  cases  at  one  time  to  find 
of  Equation  178. 


Because  the  restricted  G  is  now  not  exact,  and  even  depends  on  the  starting  values  77^  , 

it  is  necessary  to  solve  Equation  171  using  once  more  an  iterative  technique.  Noting  that  in 
1 

n  -  N  points  s  vanish,  and  applying  again  the  shorthand  Notation  180,  the  so-called  VIM 
procedure  for  a  typical  iteration  step  i  to  i  +  1  can  be  written  most  concisely  in  the  form 


,{i  +  1) 


(181) 


where  d  is  a  Boolean  diagonal  matrix  eliminating  those  characteristic  points  at  which  G  is 
applied.  For  further  details  see  Reference  24. 


c.  Example  of  Cracked  Plate 


An  interesting  example  for  the  application  of  VIM  is  a  rectangular  plate  with  a  side 
ratio  2:1,  which  is  cracked  across  the  middle  third  of  its  width  and  which  is  loaded  in  the 
longer  (axial)  direction  with  a  uniform  tensile  stress  o-Q  .  Figure  31  shows  a  quarter  of  the 
plate  and  its  idealization  with  TRIM  6  elements  We  have  used  the  cr-77  diagram  of  a 
2024- T3  aluminium  alloy,  which  can  be  closely  represented  by  the  modified  Ramberg-Osgood 
relationship, 


V  ~ 


1'10^ 
m  E 


(182) 


£  t\  O 

with  E  =  11,4  xlO  Ibf/in  ,  <^=  34,500  lbf/in  and  an  exponent  m  =  10.  Computer  applications 
confirm  that  this  problem  can  be  solved  for  the  exponential  stress-strain  Law  182  only  if  the 
VIM  technique  is  applied.  Figure  32  show,  for  various  values  of  the  tensile  loading  G  ,  typical 
curves  of  constant  normal  stress  Cyy  and  equivalent  stress  cr ,  and  also  the  plastically 
deformed  zone. 
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Figure  31.  Idealization  of  Cracked  Plate 
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Figure  32.  Plate  with  Crack 
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SECTION  V 

SOME  ASPECTS  OF  LARGE  DISPLACEMENT  ANALYSIS 

1.  GENERAL  CONSIDERATIONS 

We  presented  in  Reference  1  a  review  of  the  progress  achieved  until  autumn  1965  in  the 
application  of  the  natural  mode  technique  to  the  large  displacement  analysis  and  secondary 
instability  effects,  like  snap-through,  in  elastic  systems.  Since  then  considerable  advances 
have  been  realized  at  ISD  and  Imperial  College,  mainly  in  a  new  philosophical  direction, 
involving  continuously  varying  curvilinear  elements;  see  Section  I,  Figure  6.  The  relevant 
theory  is  very  general,  but  fits  neatly  into  the  ASKA  scheme.  Essential  refinements  of  the 
large  displacement  analysis  of  Reference  1  will  be  discussed.  These  presume  the  natural 
mode  technique  and  were  evolved  by  Kelsey  at  Imperial  College  in  1967.  Independently,  they 
were  also  developed  by  the  senior  author  in  his  lectures  at  U.C.L.A.,  and  M.I.T.  in  1966-67. 
Graduate  students  at  I.C.  working  under  the  direction  of  Kelsey  and  at  M.I.T.  under  the  senior 
author  applied  it  successfully  to  a  number  of  examples.  Some  of  the  relevant  theory  and 
illustrations  may  be  found  in  Reference  26.  The  analysis  is  essentially  based  on  a  consistent 
application  of  the  incremental  equations  (V,  24,  25)  of  Reference  1,  which  we  rewrite  in  the 
form 

R  =  K  r  =  [k  +  K  1  r  ( 1 83) 

A  A  L  E  G  J  A 

where  R^,  are  the  incremental  load  and  displacement  vectors,  K  £  the  elastic  and  K  ^  the 
so-called  geometric  stiffness  of  the  assembled  system.  For  the  sake  of  simplicity,  the 
contribution  of  the  initial  loads  is  ignored  in  Equation  183.  Given  R^(r^  ),  we  may  derive,  in 
general,  the  associated  (  R  )  and  establish  consequently  the  nonlinear  load-displacement 
relationship.  At  the  end  of  any  step  we  may  determine  first  the  natural  and  then  the  cartesian 
load  vector  P  and  P  .  The  important  point  is  now  that  Equation  183  neglects  so-called 
higher  order  terms  A  *n  equilibrium  condition  of  an  element.  This  may  lead  to 

a  considerable  deviation  from  the  true  load-displacement,  even  when  the  increments  are  small 
(but  inevitably  finite),  especially  in  the  neighbourhood  of  Instability  or  near  instability  states; 
see  Figure  33. 
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To  circumvent  this  difficulty  it  is  fairly  evident  that  we  must  apply  an  iterative 
procedure.  Consider  to  this  purpose  a  typical  loading  and  deformation  step.  The  fundamental 
transformation  matrix  4  ^  -  relating  the  natural  or  pure  deformation  vector  with  the 
total  displacement  vector  via 

P  11041 

varies  over  any  finite  increment.  We  denote  its  value  at  the  beginning  and  end  of  the  step  by 
and  respectively.  Following  Equation  (V,  10)  of  Reference  1,  the  equilibrium 

condition  would  be  exactly  satisfied,  if 


p 

9 

*  K 

pJg 

for  the  g  th  element 

(185) 

R 

II 

P  =  g  [*N  PN  5g 

for  the  complete  system 

(186) 

were  true.  Since  first  approximation  is  used  in  the  linearized  increment,  Equations  185  and 
186  hold  only  approximately.  There  remain,  in  fact  unbalanced  forces 

Ru  =  R  -  1*1  P9  (187) 

We  should  compute  a  new  linearized  increment  based  on  the  new  geometry  at  e  and  the 
residual  vector  and  continue  this  relaxation  until  Ru  is  eliminated  to  any  preset 

accuracy,  e.g.  10-3  kg.  In  Equations  185,  186  and  187  the  index  A  is  omitted  for  economy 
in  subscripts. 

Clearly  we  can  refine  and  accelerate  the  technique  by  allowing  a  second  or  even  fourth 
order  accuracy.  To  account  for  second  order  terms  we  rewrite  Equation  (V,  13)  (Reference  1) 
for  the  gth  element  in  the  evident  form 


^NgA 

2[aNb+aNe]g  ^gA 

(  188) 

^NgA  ~  *Ng  ^NgA 

~Z  kNg[aNb  +  dNelg  ^gA 

(189) 

which  means  that  the  natural  load  vector  PN^  of  the  elements  derives  from  an  average  value 
3.  N  at  the  beginning  and  end  of  the  step.  In  general.  Equations  188  and  189  prove  sufficiently 
accurate  for  a  fast  convergence  of  the  iteration  procedure  to  eliminate  the  unbalanced 
loads  Ru . 
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However,  we  may  account  for  fourth  order  terms  and  the  slope  of  the  r-R  curve  at  the 

initiation  and  end  of  the  current  step  by  remembering  Equation  (V,  19)  (Reference  1)  for  the 

definition  of  the  geometrical  stiffness  via  the  v  modal  component  matrices  A  ^  where 

o 

aa.f  * 

A  =  -~N-—  (190) 

8  dp 

We  find  for  the  higher  approximation  of  the  submatrix  g  of  (see  Equation  V,  8, 
Reference  1) 


dNgS  =  T  [aNSb  +aN8e]g  +  72  Pq  [  A8b  ”  A8e]g  0911 

and  hence 

^NgS  =  T  [aN8b  +aN8e]g/>g  +  IT  ^  [ASb  “  ASe  ]g/>g  ( l92) 

from  which  we  derive  a  more  accurate  P..  A  and  so  on. 

N  A 

hi  practical  applications  it  is,  in  general,  also  necessary  to  account  for  the  effect  of 
the  effect  of  the  PN  in  modifying  the  elastic  stiffness.  This  is  the  so-called  end  load  effect 
in  beam  elements.  Its  influence  is  described  by  the  “natural  geometrical  stiffness”  kNG  of 
Reference  1,  which  must  be  added  to  the  current  pure  elastic  stiffness  kNE  ;  see  Equation 
(V,  42)  (Reference  1).  The  theory  of  continuous  two-  and  three-dimensional  systems  using 
curvilinear  triangular  and  tetrahedron  elements,  as  developed  in  Reference  6,  includes 
ab  initio  this  effect. 

2.  EXAMPLES 

We  consider  here  two  simple  examples  for  which  exact  or  near  exact  analytical  results 
are  available.  The  first  investigates  the  large  displacements  of  a  cantilever  beam  under  a 
transverse  load  R  at  the  tip.  Figures  34  and  35  reproduce  the  nondimensional  tip  deflexions, 
Wj/t  and  -  UT  /  H,  .  They  are  based  on  a  slenderness  ratio  X  ~i/\  =120  and  presume  a 
ten-element  idealization.  True  equilibrium  was  assumed  to  be  attained  when  I0"4  kg 

hi  the  range  r£2/EI  >  1.0  the  nondimensional  loading  increment  i.2  /  El  was  fixed  at  0.1; 
for  r2-Z/E1  >1.0  the  increment  was  increased  to  0.5,  although  1,0  proves  equally  satisfactory. 
The  exact  results  are  based  on  an  infinite  direct  load  stiffness  or  \=  m .  We  observe  excellent 
agreement  to  which  contributes  the  very  small  effect  of  the  end  load  on  the  elastic  stiffness. 
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The  second  example  is  mainly  concerned  with  the  post-buckling  regime  of  a  strut  or 
beam  under  end  load;  see  Figure  36.  Here  it  is  impossible  to  achieve  a  post-buckling  equilib¬ 
rium  condition  unless  the  above  iterative  method  or  a  related  technique  is  applied.  The  effect 
of  the  natural  geometrical  stiffness  k  ^  is,  in  contrast  to  the  first  example,  significant. 
It  was  found  that  its  inclusion  yielded,  for  a  ten-element  idealization,  a  buckling  load  to 
within  an  accuracy  of  six  digits;  the  value  of  X  is  again  120.  Since  the  deflection  is  at  the 
instability  state  indeterminate,  it  is  necessary  to  impose  a  small  deflection  at  this  critical 
point.  We  may  ensure  this  by  applying  a  small  lateral  load,  which  is  removed  once  the  post- 
buckling  regime  is  secured.  Alternatively,  we  may  prescribe  for  the  beam  a  small  initial 

curvature.  The  ten  element  idealization  proves  again  timewise  economical  and  accurate.  The 

-3 

maximum  deviation  allowed  from  equilibrium  was  set  at  10  Kg;  this  measure  may  be  in¬ 
creased  without  significant  loss  of  accuracy.  For  the  non-dimensional  loading  increments 
R^JL  /  El  ,  we  chose  over  most  of  the  range  the  value  of  0.5.  However,  prior  to  reaching 
the  theoretical  buckling  load  this  was  first  reduced  to  0.1  and  then  to  0.01  and  0.005.  In  the 
post-buckling  regime  the  value  was  first  increased  to  0.01,  then  to  0.1  resp.  0.2,  the  ultimate 
value  being  once  more  0.5.  We  restrict  our  considerations  to  a  beam  with  an  initial  curvature 
of 

WQ  =  0.  001  (  I  -  cos  7T  x  /JL  )  (193) 

Figure  36  reproduces  some  of  the  results  for  the  tip  deflection  and  compares  them  with 

the  von  Mises  theory  (valid  for  small  deflections  prior  and  subsequent  to  buckling)  and  the 
Timoshenko  theory  (valid  for  large  deflections).  Excellent  agreement  is  observed,  which  also 
holds  for  high  values  of  R/R^  not  shown  in  the  diagram.  If  initial  curvature  is  varied,  we  find 
in  accordance  with  analysis  that  the  effect  is  pronounced  prior  to  and  just  beyond  the  theo¬ 
retical  buckling  state,  but  once  R/R^  Is  large  its  influence  becomes  small. 

When  an  initial  transverse  load  is  chosen,  instead  of  the  preset  curvature,  the  compu¬ 
tations  are  equally  straightforward. 

The  theory  discussed  here  in  broad  terms  has  been  extended  to  snap-through  problems 
of  arches,  spherical  caps  and  general  shells  of  revolution,  tt  is  found  that,  with  increasing 
complexity  of  the  problems,  the  new  approach  developed  in  Reference  6  proves  more  fruitful. 
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Figure  33.  Nonlinear  Behaviour  of  System- 
Softening  Effect  with  Large 
Displacements  (Due  to  Applied 
Compression  Load  P  on  Beam) 


Figure  34.  Tip  Displacements  versus 


Timoshenko 


Figure  35.  Cantilever  Beam  under  Figure  36.  Beam  with  Initial  Curvature 

Transverse  Load 
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